© ESO, 2015

1. Introduction

Our Galaxy is pervaded by an interstellar magnetic field of a few microgauss, which fills the entire disk and halo. This magnetic field manifests itself in a variety of ways, including Zeeman splitting of atomic and molecular spectral lines, Faraday rotation of polarized radio signals, synchrotron emission from relativistic electrons, and polarization of starlight and thermal dust emission. With a pressure larger than the thermal pressure of all phases and comparable to that of the cosmic rays (Cox 2005), the Galactic magnetic field (GMF) plays a crucial role in the ecosystem of our Galaxy. In conjunction with gravity, it governs the structure and the dynamics of the interstellar medium (ISM), regulates the process of star formation, accelerates cosmic rays, and channels their trajectories to confine them to the Galaxy. In addition to a large-scale regular, or coherent, component and a fluctuating component produced by interstellar turbulence (with scales up to 100 pc; e.g., Gaensler & Johnston 1995; Haverkorn et al. 2008), the GMF also possesses an ordered random (e.g., Beck 2009; Jaffe et al. 2010), or striated random (Jansson & Farrar 2012a), component, whose orientation remains nearly constant over large scales, but whose strength and sign vary on small scales. Such fields are probably produced through compression or shearing of isotropic random fields by the Galactic differential rotation, or at large-scale spiral arm shocks, or else by rising hot plasma bubbles.

Our knowledge and understanding of the GMF has improved considerably over the past few years, as a result of both progress in the quality (sensitivity and resolution) of radio observations and extensive modelling efforts (e.g., Sun et al. 2008; Sun & Reich 2010; Ruiz-Granados et al. 2010; Jaffe et al. 2010, 2011; Pshirkov et al. 2011; Fauvet et al. 2012, 2013; Jansson & Farrar 2012a,b). However, the existing radio observations have inherent limitations, as both Faraday rotation measures (RMs) and synchrotron (total and polarized) intensities are quantities integrated over the line of sight (LOS), which depend on the poorly constrained density distributions of thermal and relativistic electrons, respectively. A promising avenue to obtain a more complete and more robust picture of the GMF structure is to complement the radio data with Planck¹ measurements of the polarized thermal emission from interstellar dust, which is independent of the electron densities.

A glance at the Planck all-sky intensity maps (Planck Collaboration I 2014) reveals that, in addition to the mottled structure of the cosmic microwave background (CMB) at high Galactic latitudes, the dominant pattern is that of the emission from our Galaxy. At the lowest frequencies, from the 30 GHz to 70 GHz bands of the Planck Low Frequency Instrument (LFI, Bersanelli et al. 2010), synchrotron emission dominates; at the highest frequencies, from the 217 GHz to 857 GHz bands of the High Frequency Instrument (HFI, Lamarre et al. 2010), thermal emission from interstellar dust is the dominant emission. These foregrounds have to be understood and taken into account for detailed CMB studies, but they also provide a unique opportunity to study the Galaxy’s ISM.

In particular, the thermal dust emission is linearly polarized (e.g., Benoît et al. 2004; Vaillancourt 2007). This polarized emission overpowers any other polarized signal at the higher Planck frequencies (e.g., Tucci et al. 2005; Dunkley et al. 2009; Fraisse et al. 2009). In addition to hindering the detection of the sought-after, odd-parity, B-mode polarization of the CMB (Planck Collaboration Int. XXX 2015), the polarized dust emission provides, in combination with the emission spectrum itself, a powerful constraint on the physical properties of the dust and on the structure of the magnetic field in the Galaxy.

The linear polarization of the thermal dust emission arises from a combination of two main factors. Firstly, a fraction of the dust grain population is non-spherical, and this gives rise to different emissivities for radiation with the electric vector parallel or orthogonal to a grain’s longest axis. Secondly, the grains are aligned by the interstellar magnetic field because they are rotating, probably with differing efficiencies depending on grain size and composition (Draine & Fraisse 2009). While the details of this process remain unclear (Lazarian 2003, 2007), there is a consensus that the angular momentum of a grain spun up by photon-grain interactions (Dolginov & Mitrofanov 1976; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007; Hoang & Lazarian 2008) becomes aligned with the grain’s shortest axis, and then with the magnetic field via precession (e.g., Martin 1971). The end result is that, if we look across magnetic field lines, the rotating grain has its long axis orthogonal to the field lines, and accordingly dust emission is linearly polarized with its electric vector normal to the sky-projected magnetic field².

A related phenomenon occurs at near-UV/optical/NIR wavelengths (e.g., Martin 2007), where the light from background sources becomes linearly polarized as a result of dichroic extinction by the aligned dust grains (Davis & Greenstein 1951). Because extinction is higher for light vibrating parallel to the grain’s longest axis, i.e., perpendicular to the field lines, the transmitted light is linearly polarized with its electric vector parallel to the sky-projected magnetic field. In fact, historically, the optical polarization caused by dust extinction led to the prediction that thermal dust emission would be polarized in the millimetre and submillimetre domains (Stein 1966). The predicted orthogonality of the electric vectors in the optical and submillimetre on the same line of sight has been demonstrated (Planck Collaboration Int. XXI 2015).

Thus, polarized thermal dust emission carries important information on the interstellar magnetic field structure, on the grain alignment mechanisms, and on the grain geometrical and physical properties. For example, polarization observations between 300 μm and 3 mm, essentially the domain of the Planck HFI instrument, can potentially discriminate between the polarizing grain materials, e.g., silicate and graphite dust versus silicate-only grains (Martin 2007; Draine & Fraisse 2009; Planck Collaboration Int. XXI 2015; Planck Collaboration Int. XXII 2015).

The far-IR dust thermal emission being a tracer of the dust mass along the LOS, sensitivity limits explain why detailed dust polarized emission was observed mostly in fairly dense, massive regions of the ISM (Dotson et al. 2000; Curran & Chrysostomou 2007; Matthews et al. 2009; Dotson et al. 2010), in general close to the Galactic plane. Measurements of the more diffuse medium were obtained at relatively low (2deg) angular resolution. At these large scales, the Archeopsballoon experiment (Benoît et al. 2004; Ponthieu et al. 2005) detected the thermal dust emission polarization at 353 GHz. The highest frequency channel of WMAP (Page et al. 2007; Bennett et al. 2013), 94 GHz, picked up the long-wavelength Rayleigh-Jeans tail of the diffuse dust emission and its polarization (in addition to synchrotron emission).

The Planck satellite’s HFI instrument has led to the first all-sky survey of the polarized submillimetre and millimetre sky, where thermal dust emission dominates. At 353 GHz, the Planck data have an angular resolution of 5′. The polarization sensitivity was expected to be such that, at a resolution of 15′, ISM structures with A_V = 1 mag would be detected with a relative uncertainty on the polarization fraction of about 40% and an uncertainty on the polarization angle of about 30° (Pelkonen et al. 2009). These figures improve significantly at higher A_V and/or lower resolution. The polarized Planck data bring the first all-sky fully sampled map of the polarized emission from dust. As such, they provide unprecedented information on the magnetic field geometry and the dust polarization properties relevant to the disk of the Milky Way (MW) and star forming regions, for which they provide statistical information that is missing in stellar polarization extinction data. It should be emphasized, however, that the dust polarized emission provides information mostly on the orientation of the sky-projected magnetic field and only very indirect indication about the angle of that field with respect to the plane of the sky, and it is expected to be insensitive to the field strength.

This paper presents a subset of the Planck polarization data and their large-scale statistical properties. A companion paper (Planck Collaboration Int. XX 2015) analyses the variations of the polarization fraction and angle described here, in the framework of simulations of anisotropic magneto-hydrodynamic (MHD) turbulence. Two other papers provide a detailed analysis of the wavelength dependence of the dust polarization, as seen by the HFI instrument (Planck Collaboration Int. XXII 2015) and a comparison between the dust polarization at visible and submillimetre wavelengths (Planck Collaboration Int. XXI 2015).

In Sect. 2 we describe the data, including discussion of systematic effects and the effects of the CMB intensity and polarization. Maps are presented in Sect. 3, as well as the statistics of the data. Sect. 4 discusses the implications of the 353 GHz polarimetry for our understanding of the GMF structure, and the conclusions are drawn in Sect. 5. Three appendices discuss the smoothing of the noise covariance matrices, which is needed when the original data are averaged, the debiasing methods for obtaining polarization estimates, and tests for the effects of systematic noise bias on the structures that we observe in maps of the polarization angle dispersion function.

2. Data

The Planck mission results are presented in Planck Collaboration I (2014) and the in-flight performance of the two focal plane instruments, the High Frequency Instrument (HFI) and the Low Frequency Instrument (LFI), are given in Planck HFI Core Team (2011) and Mennella et al. (2011), respectively. The data processing and calibration of the HFI data used here are described in Planck Collaboration VI (2014), Planck Collaboration VII (2014), Planck Collaboration VIII (2014), Planck Collaboration IX (2014) and Planck Collaboration X (2014). The data processing and calibration of the LFI data are described in Planck Collaboration II (2014), Planck Collaboration III (2014), Planck Collaboration IV (2014), and Planck Collaboration V (2014).

The Planck polarization and total 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 solely from total intensity with polarization data at small scale (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, 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. The error-bars we quote include uncertainties associated with residual systematics as estimated by repeating the analysis on different subsets of the data. We have also checked our data analysis on the latest version of the maps available to the collaboration, to check that the results we find are consistent within the error-bars quoted in this paper.

The maps used include data from five independent consecutive sky surveys (called Survey1-Survey5) for HFI, taken six months apart. Due to the scanning strategy of the Planck mission, surveys taken one year apart (i.e., odd surveys 1 and 3 and even surveys 2 and 4) share the same observing pattern, which is different for even and odd surveys. Survey5 had a different scan pattern from the other odd-numbered surveys, owing to a change in the precession phase. The products also include data binned into the first and second halves of the Planck stable pointing periods, or “half-rings” (called HR1 and HR2). Both single-survey and half-ring data are used for consistency checks and to assess the level of systematic effects. Here, we only analyse the polarization data at 353 GHz, which is the highest frequency Planck channel with polarization capabilities and the one with the best S/N for dust polarization. We use the 30 GHz LFI data in our comparison of the dust emission at 353 GHz with the microwave and radio synchrotron emission presented in Sect. 4.4.

In the Planck map-making process (Planck Collaboration VIII 2014), measurements from various detectors at the same frequency are combined to obtain the Stokes parameters (I, Q, and U) at each position on the sky. The reconstructed polarization is a linear combination of the weighted differences between the signal from pairs of polarization sensitive bolometers (PSBs) with different orientations on the sky. The resulting maps of the Planck Stokes parameters Q and U used in this paper are shown in Fig. 1. The corresponding map of the observed polarization intensity P = (Q² + U²)^{1 / 2} is shown in Fig. 2. The total intensity map used in this work is shown in Fig. 5.

2.1. Conventions and notations

The relations between the observed Stokes parameters (I, Q, and U) and the polarization fraction (p) and polarization angle (ψ) are given by (1)and (2)where the two arguments function arctan(Y,X) is used to compute atan(Y/X) avoiding the π ambiguity, such that (3)For the Stokes parameters provided in the Planck data, the angle convention above is with respect to Galactic coordinates with −90° <ψ < + 90°, ψ = 0° toward Galactic north, and ψ becoming positive toward Galactic west, the direction of decreasing Galactic longitude (i.e., ψ increases clockwise). Note that this convention is the one used in the HEALPix³ software (Górski et al. 2005), but is different from the IAU convention (Hamaker & Bregman 1996), which is ψ = 0° toward Galactic north but with ψ becoming positive toward Galactic east, the direction of increasing Galactic longitude (i.e., ψ increases counterclockwise). The conversion between Planck Stokes parameters and the IAU convention is given by: (4)In this paper, all quoted values of the polarization angle are given in the IAU convention.

2.2. Bandpass mismatch leakage correction

Owing to the way the polarization maps are constructed, any instrumental difference between detectors of the same channel may produce a fake polarization signal, even for unpolarized sky signal inputs. This is the case for the bandpass mismatch (BPM) between detectors that affects Planck polarization maps. In practice, the effect corresponds to a leakage term from total intensity I into polarization Q and U. The BPM polarization leakage effect is therefore strongest in regions of high intensity, i.e., along the Galactic plane, and affects both p and ψ. Because the 353 GHz intensity data used here are calibrated on the CMB signal, no BPM leakage is produced by the CMB anisotropies. Other astrophysical emission sources, however, produce BPM polarization leakage.

Knowing the actual Planck sky scanning strategy and the orientations of the polarization sensitive bolometers in the focal plane, the BPM polarization leakage corrections can be estimated from the relative responses of each detector to a given sky astrophysical emission. The Planck Collaboration is exploring different methods to compute the relative responses of detectors, as well as to produce total intensity maps for each sky emission source. Two methods have been used to determine the relative responses (Planck Collaboration IX 2014). The first one (method A) involves computing the BPM leakage between bolometers using the ground-measured bandpasses (Planck Collaboration IX 2014). The second one (method B) deduces the relative detector response on regions of the sky where we can obtain I, Q, and U maps for each detector individually. Note that this can only be performed in limited regions of the sky, outside the Galactic plane, which have been scanned in a large number of configurations, allowing for the full reconstruction of I, Q, and U per detector. A comparison between the two methods is presented in Planck Collaboration IX (2014).

When folding the above coefficients into the Planck scanning strategy, we have chosen to produce template maps of the BPM leakage contribution for each frequency (ν) channel, for each bolometer (b(ν)) and for each Stokes parameter (X being Q or U). The BPM polarization leakage correction is (5)where R_b(ν) represents the detector relative responses and I_ν is the sky total intensity. For the purpose of the study presented here we only take into account BPM leakage from dust thermal emission, because this is the dominant term at 353 GHz. The template maps in Eq. (5) were computed using the Planck thermal dust model described in Planck Collaboration XI (2014). We used the standard Planck map-making procedure presented in Planck Collaboration VIII (2014). Note that the Planck 353 GHz channel also includes emission from the CO (J = 3 → 2) line (see Planck Collaboration VI 2014), which should also in principle be included in the BPM leakage correction. This, however, is relatively weak with respect to dust thermal emission and the corresponding BPM effect is expected to be small compared to that from dust. Because we do not concentrate on regions with strong molecular emission in this paper, no correction was applied for the CO emission BPM leakage.

Figure 3 shows the effect of the correction for BPM on the observed distribution of polarization angles toward the plane of the Milky Way (| b_II | < 5°) in the four Galactic quadrants (Q1, Q2, Q3 and Q4, defined by 0° <ℓ_II < 90°, 90° <ℓ_II < 180°, 180° <ℓ_II < 270°, and 270° <ℓ_II < 360°, respectively). When no BPM leakage correction is applied, angles are observed to be distributed around + 20° and −5° for the inner (Q1 and Q4) and outer (Q2 and Q3) MW regions, respectively. The difference in sign is due to the difference in average detector orientation during Galaxy crossings, resulting from the relative orientation of the scanning strategy and the Galactic plane. Using the two methods discussed above for the determination of the coupling coefficients leads to similar BPM leakage estimates. Note also that because the magnetic field is expected to be statistically aligned with the Galactic plane (see, e.g., Ferrière 2011) we expect the polarization direction toward the plane to be on average around ψ = 0°. The fact that both correction methods bring the peak of the histograms toward this value confirms the validity of the BPM correction method used here. In the following, we adopted the coefficients from method B. We note, however, that although the situation is clearly improved by the BPM leakage correction, the average observed angle distributions still peak a few degrees away from ψ = 0°, with the same sign pattern as for the uncorrected data. This could in principle be due to incomplete correction. However, preliminary tests have shown that the remaining correction could be due to non-linearity in the analogue-to-digital conversion (ADC) of the signal, which produces an additional correction with the same sign as observed here and roughly the right amplitude.

We do not attempt here to fully assess the quality of the different corrections, but simply use them to estimate where on the sky the uncertainties in the corrections are small enough to be unimportant for this study. A plot of the BPM-leakage-corrected polarization angle versus the uncorrected polarization angle shows the magnitude of the correction, while the correlation coefficient gives a quantitative measure. For the different corrections considered above, the correlation coefficient is over 0.95 for most regions of the sky at | b_II | > 5°. Above | b_II | = 10°, the correlation coefficients are above 0.98, implying that the correction becomes very small. This is a natural result of the fact that the intensity that is leaking into polarization is brightest toward the Galactic plane. As measured from the difference between method A and B, the corresponding uncertainties on the polarization angle ψ and fraction p are | Δψ | < 10° and Δp < 1%, respectively, toward the inner Galactic plane. These uncertainties become less than the random errors away from the plane. However, BPM leakage corrections are probably not the dominant uncertainty at high Galactic latitudes and very low signal levels, where other systematic effects remaining in the data become more important (see Sect. 2.4). For this reason, we do not discuss specifically the polarization properties in the lowest brightness sky area in this paper and defer this discussion to future papers.

The above discussion applies to the HFI data, but we will also compare the thermal dust emission at 353 GHz to the 30 GHz emission from LFI, which has a similar bandpass leakage issue. The LFI BPM correction is discussed in Planck Collaboration II (2014), where the principle difference is the presence of multiple astrophysical foregrounds, with different spatial and spectral distributions. The component separation products are therefore used in the LFI BPM correction. From a comparison of the different surveys, we estimate that the uncertainties are of the order 10 μK in the polarized intensity and dominated by the noise rather than the leakage except in the innermost plane (| ℓ_II | < 30° and | b_II | < 3°), where the effect is only slightly above the noise level. For the polarization angle, we estimate the uncertainties as roughly 15° in the plane (| b_II | < 5°) and 35° away. Again the uncertainty appears dominated by noise, with no obvious structure related to the bandpass leakage or scan pattern. We have also cross-checked with WMAP 23  GHz data and verified that the results in Sect. 4.4 are very similar.

2.3. Deriving polarization parameters

The polarization parameters I, p, and ψ are derived from the observed Stokes parameters I, Q, and U using the Bayesian method described in Montier et al. (2015a). This method extends that described in Quinn (2012) by using the full 3 × 3 noise covariance matrix of each pixel. The effective ellipticity, as defined in Montier et al. (2015a), characterizes the shape of the noise covariance matrix and couples all the terms in Q and U. ϵ_eff = 1 corresponds to the case described in Quinn (2012), whereas ϵ_eff> 1 means that the relation between C_QQ,C_QU,C_UU is not trivial, and there are asymmetries in the noise covariance matrix. We calculated ϵ_eff for the Planck data used here. At 1° resolution it is normally distributed with a mean value of 1.12 and a standard deviation of 0.04. At the full Planck resolution, the distribution of ϵ_eff is a bit wider (standard deviation of 0.05), but the mean value does not change. Thus, although they are not very strong, the asymmetries of the noise covariance matrix cannot be neglected, and the Bayesian method is well suited for the analysis of the data.

We use a flat prior on all three parameters p, ψ, and I over a range centred on the conventional value given by Eqs. (1) and (2) for p and ψ and the observed value for I, and a width corresponding to 20σ, where σ is the conventional estimate for the uncertainties (see Appendix B.1). The range on p and ψ is further limited to −1 <p < 1 and −90° <ψ < 90°, respectively, allowing negative values of p in order to reduce bias in the posterior probability. We compute the 3D posterior probability distribution function (PDF) using 2⁷ values on each axis over the parameter range. The values of the polarization parameters are obtained using the mean posterior (MP) estimator on the posterior 3D PDF. A comparison between the polarization parameters and uncertainties obtained with this method and using the conventional approach described in Appendix B.1 is shown in Fig. B.1 for the Planck data at 1° resolution.

When spatial smoothing is applied to the polarization data, Stokes parameter maps are convolved with a Gaussian kernel of the appropriate width using the dedicated smoothing software part of the HEALPix library, which guarantees proper transport of Q and U. The maps are then resampled to larger pixel size (as specified by the HEALPix N_side parameter) so as to preserve full sampling of the data (pixel size smaller than 1/2.4 times the data FWHM resolution). The corresponding smoothing of data covariances was performed using the method described in Appendix A. The corresponding smoothed maps of p and ψ are then computed as described above. The statistical uncertainties in p and ψ ( and , respectively) have been estimated as described in Appendix B.3.

2.4. Impact of systematic effects, CIB, ZL, and CMB

We assessed the level of contamination by systematic effects by comparing the maps of p and ψ obtained at 1° resolution for the full Planck data with those obtained for the various individual Planck surveys (see introduction to Sect. 2). We constructed maps of systematic uncertainties on p and ψ ( and , respectively) by averaging these differences over the Planck individual surveys. These were added to the statistical uncertainty maps and , to obtain the total uncertainty maps used in the rest of the analysis.

In this paper we only show the Planck polarization data and derived quantities where the systematic uncertainties are small and where the dust signal dominates total emission. For this purpose we defined a mask such that and I₃₅₃> 0.1 MJy sr^-1. We defined the mask at a resolution of 1° and smoothed it to 3° resolution to avoid complex edges. As a result, the maps shown exclude 21% of the sky. Note that a different mask is used for the polarization angle dispersion function, as defined in Sect. 3.3.

The cosmic infrared background (CIB) is due to emission from a large number of distant galaxies with random orientations and is expected to be, on average, unpolarized. However, it can contribute non-negligible emission at 353 GHz in low brightness regions of the sky and hence reduces the apparent degree of dust polarization. The zero level of the 353 GHz total intensity map has been established by correlation with Galactic H i, using the method described in Planck Collaboration XI (2014), as was done for the publicly released 2013 maps. This offset is 0.0887 MJy sr^-1 (uncertainty 0.0068 MJy sr^-1) and was subtracted from the total intensity map we use, which therefore does not contain the CIB monopole. We added the corresponding uncertainty in quadrature with the uncertainty of the total intensity, so that the statistical uncertainties on p include the uncertainty on the CIB subtraction.

The zodiacal light (ZL) has a smooth distribution on the sky. From the model constrained by its detection in the Planck bands (Planck Collaboration XIV 2014), its median total intensity at 353 GHz is 1.9 × 10^-2 MJy sr^-1 over the sky area studied here, and reaches ≃4.3 × 10^-2 MJy sr^-1 in dust lanes near the ecliptic plane. Its polarization in the submillimetre is currently unconstrained observationally. Because this intensity is subdominant over most of the sky fraction and the polarization level of ZL is currently unknown, we apply no correction for the possible contribution of ZL. We note that, if ZL was assumed unpolarized, subtracting its intensity would raise the observed polarization levels by about 0.5% of the observed polarization fraction, on average over the sky region studied here, and would not change the observed polarization angles. We have checked that no noticeable systematic variation of the polarization fraction is detected in our maps along zodiacal dust lanes.

CMB fluctuations are polarized at a level of 0.56 mK (Kovac et al. 2002) at a resolution of about 1deg, which corresponds to 1.6 × 10^-4 MJy sr^-1 at 353 GHz. In the mask we use here, the effect of CMB polarized fluctuations is therefore negligible and we did not attempt to correct for those fluctuations.

No additional correction was applied to the data.

2.5. External data

In Sect. 4.4, we compare the Planck HFI polarization maps with low-frequency radio and microwave observations that are dominated by synchrotron emission over most of the sky. These include:

• the 408 MHz total intensity map of Haslam et al. (1982) from the LAMBDA⁴ site;

• the 1.4 GHz total intensity map of the northern (Reich 1982; Reich & Reich 1986) and southern (Reich et al. 2001) sky;

• the 1.4 GHz polarized intensity maps of the northern (Wolleben et al. 2006) and southern (Testori et al. 2008) sky.

For the analysis in Sect. 4.4, the Planck HFI and LFI maps are smoothed to 1° FWHM resolution to match these radio data and downgraded to N_side = 256. Most of the 1.4 GHz maps are available on the Bonn survey site⁵ as FITS images in Cartesian coordinates. They are converted into HEALPix using the procedure described in Paradis et al. (2012) and are made available in this form on the CADE site⁶. The resolution of the observations is roughly 1°, and so no additional smoothing is applied to the radio data. The total intensity map at 1.4 GHz is estimated to have an offset of 2.8 K (Reich et al. 2004) due to the combination of zero-level calibration uncertainty, unresolved extragalactic sources, and the CMB, and so this was subtracted from the data.

The total intensity data include thermal bremsstrahlung (free-free) emission, particularly in the plane. This is not negligible at 408 MHz or 1.4 GHz. We use the WMAP MEM free-free solution (Gold et al. 2011) to subtract it. We note that this free-free template likely includes anomalous dust emission, and there are indications that it is an overestimate by roughly 20 to 30% (Alves et al. 2010; Jaffe et al. 2011). Because synchrotron dominates over free-free emission at low radio frequencies, even on the Galactic plane, the uncertainties on the free-free correction are not expected to affect the qualitative comparison with dust emission in this paper. But the MEM template is not sufficiently accurate to correct for free-free when the synchrotron is subdominant at 30 GHz. Furthermore, the 30 GHz total intensity also includes anomalous dust emission for which we have no correction. We therefore do not use 30 GHz in total intensity, but only in polarization.

3. Description of the Planck polarization maps

Figure 4 shows the maps of the polarization fraction (p) at a resolution of 1°. The map in Fig. 5 is based on the polarization direction, also at a resolution of 1°. Both figures also show the corresponding map of the total uncertainty, which includes the contribution from statistical and systematic uncertainty estimates, as described in Sect. 2.4. The maps were masked as described in Sect. 2.4 in regions where large residual systematic uncertainties were evident or where the total intensity at 353 GHz is not dominated by dust emission. Figures 4 and 5 were constructed using the Bayesian method described in Sect. 2.3, Montier et al. (2015a), and Appendix B.3, in particular the Mean Posterior Bayesian estimator defined in Montier et al. (2015b). These figures are discussed in Sects. 3.1 and 3.2. In Fig. 6 we highlight several regions of interest that we will discuss below; parameters of these regions are given in Table 1.

3.1. Polarization fraction

As seen from Fig. 4, the measured polarization fraction shows significant variations on the sky. One of the aims of this paper is to characterize those variations as a step toward understanding their origin. These characteristics are compared to those of polarized emission maps computed in simulations of anisotropic MHD turbulence in a companion paper (Planck Collaboration Int. XX 2015).

Figure 4 shows that the polarization fraction of the thermal dust emission can reach up to about 20% in several large-scale regions of the sky. This is particularly the case in the second Galactic quadrant (ℓ_II ≃ 145°, b_II ≃ 0°, including a region at low latitude known as “the Fan”⁷), the Perseus area (ℓ_II ≃ 143°, b_II ≃ −25°), the Loop I area (ℓ_II ≃ 40°, b_II ≃ + 45°) and a region in Microscopium (ℓ_II ≃ 15°, b_II ≃ −40°). The large-scale distribution of these regions is consistent with predictions from the Galactic magnetic field model used in the Planck Sky Model (PSM) (Delabrouille et al. 2013). This model, based on a simple description of the spiral magnetic field structure of the Milky Way, was optimized to match the WMAP and Archeopsdata (e.g., Fauvet et al. 2011, 2012). The model computes a dust polarization geometrical efficiency factor g_d, which results from the LOS integration of the magnetic field direction, weighted by the assumed dust density distribution. This factorhas a maximum toward the Galactic anticentre at a position corresponding roughly to that of the Fan region and shows two strong maxima at mid-latitude toward the inner Galaxy ℓ_II ≃ 0° and | b_II | ≃ 45° which match fairly well with the highly polarized regions detected with Planck around Aquila Rift, Libra, and Pavo, above and below the Galactic plane, respectively.

Figure 7 shows the distribution of p as a function of the polarization geometrical efficiency factor g_d, in particular the polarization fraction computed for the bi-symmetrical spiral model of the magnetic field on large scales from Miville-Deschênes et al. (2008). It can be seen that the upper envelope of the distribution roughly matches that allowed by the model, indicating that it is only when g_d is close to unity, i.e., toward regions where the large-scale Galactic magnetic field is preferentially orthogonal to the LOS, that high p values are observed. The dispersion of the points below this line is explained by small scale variations of p of a different origin, described later in this paper.

Figure 8 displays the histogram of the polarized fraction p over the sky fraction shown in Fig. 4, the whole Galactic plane (| b_II | < 5°) and the inner Galactic plane (| b_II | < 5°, | ℓ_II | < 90°) at a resolution of 1°. In the plane, the most likely value of p is a few percent while the rest of the mid-latitude sky has a wider distribution, with a peak of the histogram near 4%. The maximum p values can reach about 20%. A more accurate determination of the maximum p value p_max, taking into account the effects of data resolution and noise, is given in Sect. 4.1 and leads to a similarly high value. We note that this maximum value is much higher than values reported previously from ground-based observations in the submillimetre. This is mainly because such low brightness regions are too faint to be observed from the ground, and because higher column density and brighter regions, which can be observed from the ground, have a tendency to be less polarized than faint regions (see Sect. 4.2). However, we show in Sect. 4.2 that the high maximum p values in low brightness regions is not a trivial consequence of the decrease of the emission intensity. We also note that the high polarization fractions observed here are consistent with the value inferred from the Archeopsmeasurements at 353 GHz, which was derived to be as high as 10–20% (Benoît et al. 2004) along the outer Galactic plane, a region which includes the Fan region.

Figures 9–11 display maps around some of the regions outlined in Fig. 6 and listed in Table 1. Figure 9 shows the Aquila Rift and Fan regions, which show high polarization fraction. These highly polarized regions are generally located in low intensity parts of the sky (e.g., Microscopium, Libra, Pavo or Ara), or at the edge of bright regions (e.g., the Aquila Rift). They are also located in regions of the sky where the polarization direction is uniform over very large areas. For instance, in the Fan region, the magnetic field is oriented almost completely parallel to the Galactic plane (i.e., polarization is orthogonal to the plane) with high accuracy over a region spanning more than 30°, where the polarization fraction consistently has p> 8% and reaches p ≃ 15% in some areas. Similarly, the highly polarized Aquila Rift region has a B-field sky projection aligned with the elongated structure of the ridge and the nearby Loop I over most of the extent of the source, and the polarization fraction there reaches up to 20%. The highly polarized region is in fact located on the gradient of the dust emission of the Aquila Rift, and mid-way between the Aquila Rift itself and the radio emission of Loop I. In the Perseus region, the large polarization also appears in fairly low brightness regions, where the orientation of the field is coherent over regions of the sky with typical sizes of a few degrees. Some of these structures were detected in polarized light at other wavelengths. For instance, the Fan, Perseus, and Loop I regions seem to have counterparts detected in polarized thermal dust and synchrotron emission, as well as in the WMAP foreground emission (Gold et al. 2011; Ruiz-Granados et al. 2010; Jansson & Farrar 2012a, and references therein) and in Faraday depth surveys of polarized emission at radio frequencies, such as the Global Magneto-Ionic Medium Survey (GMIMS; Wolleben et al. 2010b). In particular, from the Faraday depth data of GMIMS, a significant portion (about 5%) of the sky was found to be dominated by the magnetic field around a nearby H i bubble (at a distance of ≃100 pc) in the general direction of the Loop I region described above (Wolleben et al. 2010a). In general, such regions are identified with nearby Galactic structures (e.g., supernova remnants and bubbles), which can distort the underlying more regular large-scale pattern of the Galactic magnetic field. Finally, other regions, such as Microscopium, have almost no known counterpart structure at other wavelengths. The area around Ara and Triangulum has been identified only as a region with warmer dust in Planck Collaboration XIX (2011). Here too, the polarization fraction is typically p> 10% (see also Sect. 4).

As seen in Figs. 4 and 8, the inner Galactic plane shows much lower polarization fractions than the highly polarized regions described above. Because the polarized emission also increases in that region, the decrease of p is a real depolarization effect. It results from a combination of depolarization due to LOS integration and the presence of dense clouds, which tend to have lower polarization fractions (see Sect. 4.2).

Note that the map of polarized intensity itself, Fig. 2, exhibits narrow features where the polarized intensity drops. These regions are also seen as narrow features where the polarization fraction drops (Fig. 6). For instance, one can be followed over several tens of degrees, rising from the Galactic plane at ℓ_II = 90deg and crossing the Polaris Flare region in Figs. 6 and 17. Inspection of Figs. 9–11 shows that these regions are not trivially peaks of the total intensity with no polarized intensity counterpart. They are sometimes regions with higher gas column density N_H (see Sect. 4.2), but not always. They can also be regions where the orientation of the field changes more abruptly (see Sect. 3.2 for a full discussion).

3.2. Polarization angle

Figure 5 shows the large-scale distribution of the polarization direction. In the figure, the direction shown by the normalized segments is that of the observed polarization direction (ψ) rotated by 90°. The figure therefore shows the orientation of the apparent magnetic field (⟨ B_⊥ ⟩). In the simplified case that the direction of B remains homogenous along the LOS, ⟨ B_⊥ ⟩ measures the projection of B onto the plane of the sky, i.e., perpendicular to the LOS. However, in the more realistic case of a disordered B structure and inhomogeneous dust distribution along the LOS, it is important to remember that ⟨ B_⊥ ⟩ is a LOS-averaged quantity, weighted by dust emission.

Figure 5 shows that, toward the Galactic plane, ⟨ B_⊥ ⟩ is mostly oriented along the plane, corresponding to a polarization angle close to 0°. This is especially the case toward the outer MW regions. There are a few exceptions, in particular toward lines of sight that are tangent to spiral arms (Cygnus X, ℓ_II ≃ 81°, b_II ≃ 0°; Carina, ℓ_II ≃ 277°, b_II ≃ −9°), where the polarization signal is actually the smallest in the plane because in those regions the large-scale magnetic field is parallel to the LOS. This was already noted by Heiles (1996) (and references therein) and Benoît et al. (2004). We also note that the homogeneity of the field orientation being parallel to the plane extends away from the plane and up to | b_II | ≃ 10° in many regions (in particular the Fan). At intermediate latitudes, the field orientation follows a few of the well known filamentary intensity structures of the local ISM. In particular, this is the case for the Aquila Rift and most of Loop I (outside the latitude range b_II ≃ 50°−60°), where the structure of ⟨ B_⊥ ⟩ follows the intensity flare and loop elongation. As addressed earlier, this orientation of ⟨ B_⊥ ⟩ in those regions was already noted in the synchrotron polarized maps of WMAP (Gold et al. 2011). Other regions, however, show a variety of relative orientations between the field projection and intensity structures, which can also be orthogonal in some instances. Thus studies with Planck submillimetre polarization (Planck Collaboration Int. XXXII 2015; Planck Collaboration Int. XXXIII 2015) hold promise as a valuable complement to optical and near infrared polarization studies of the relative orientation of the magnetic field and structure (e.g., Goodman et al. 1990; Chapman et al. 2011).

3.3. Polarization angle dispersion function

In order to quantify the regularity of the B field revealed by the polarization measurements, we use the polarization “angle dispersion function” given by (6)where Δψ_xi = ψ(x) − ψ(x + δ_i) is the angle difference between ψ(x), the polarization angle at a given sky position x (the central pixel), and ψ(x + δ_i) the polarization angle at a sky position displaced from the centre by the displacement vector δ_i. The average in Eq. (6) is taken over an annulus around the central pixel of radius δ = | δ | (the “lag”) and width Δδ and containing N pixels. Note that the angle dispersion function as defined in Eq. (6) is a two-point function, but depends both on position and lag. In that sense, it is distinct from structure functions often used to describe polarization direction (see Serkowski 1958; Kobulnicky et al. 1994; Hildebrand et al. 2009) which do not depend on position. Note also that different ways of characterizing the field geometry have been used in the literature, such as the normalized spatial gradient of the polarized intensity (| ∇P | /P) discussed in Burkhart et al. (2012) and used in Iacobelli et al. (2014). We have also estimated this quantity using the Planck data and the map obtained exhibited the same structure as the map of .

In practice, Δψ_xi is computed from the Stokes parameters as (7)where indices x and i stand for the central and displaced values, respectively. We use Δδ = δ so that individual pixels are not counted twice when estimating at different lags.

The polarization angle dispersion function measures the inhomogeneity of the polarization direction orientation, irrespective of absolute direction. It provides important information on the magnetic field distribution and orientation (see, e.g., Falceta-Gonçalves et al. 2008; Poidevin et al. 2013). Regions where the sky projection of the magnetic field is well ordered have , while regions with a twisted or disordered B field can in principle have up to . In addition, because the Planck convention for Q and U is defined with respect to the Galactic coordinate system, even a homogeneous field would produce , due to the coordinate singularity at the poles. In order to avoid this, we have rotated Q and U locally to a different coordinate system so that the current point is at the equator of the new system, before applying Eq. (7). When the signal is dominated by noise, converges to (≈ 52°). The angle dispersion function is observed to increase with δ, as the coherence is gradually lost when going further away from a given point of the sky. It is expected to increase linearly with lag in the presence of an ordered magnetic field and to steepen at small lags due to either the turbulent component of the magnetic field or the angular resolution of the data used (see, e.g., Hildebrand et al. 2009). The dependence of on lag δ can be better probed from the analysis of individual regions at higher resolution, either in emission or in absorption toward stellar clusters (Magalhães et al. 2005; Falceta-Gonçalves et al. 2008; Franco et al. 2010).

Like other quadratic functions, is biased positively when noise is present. As described in Hildebrand et al. (2009), can be debiased using (8)where is the variance on . In the conventional approach, can be expressed as a function of σ_ψ through partial derivatives as (9)where σ_ψ is the straightforward standard deviation of ψ, i.e., computed without resolving polarization angle ambiguities. However, this approximation is valid only close to the solution and leads to a poor estimate of the bias at low S/N. Nonetheless, it is clear from Eqs. (9) and (B.7) that regions with low polarization having higher values of σ_p/p have higher σ_ψ and therefore more biased .

In order to assess the importance of the bias, we use the two independent half-ring maps H₁ and H₂ of the Planck data to compute an unbiased estimate of as (10)where is the angle difference for half-ring map H_j, i.e., . In practice, is computed as (11)Although is unbiased, it suffers from higher noise because only half of the Planck data are used. Note also that, unlike , can be negative.

We evaluate from the full Planck survey (we call this estimate simply by default) and at each pixel of the map using Eqs. (7) and (11), respectively. We also perform a Monte Carlo noise simulation on I, Q, and U for each pixel using the full covariance matrix (using Eq. (A.23)), and assuming that different pixels have independent noise and that the half-ring maps have independent noise. This simulation is used to construct the PDF of , and using 1000 noise samples. We then compute the mean posterior estimates and uncertainties of , and by integrating over the PDF.

Figure 12 shows the sky distribution of computed from the full survey at 353 GHz at 1° resolution for a lag of δ = 30′. Figure 13 shows the same maps obtained from the half-ring survey correlation (), as well as for individual Planck surveys. The mask used in the these figures was obtained from the uncertainty on , , derived from the Monte Carlo analysis described above. The mask is such that the S/N on is larger than 3 () and retains 52% of the sky at the adopted analysis resolution of 1°. The differences between individual panels of Fig. 13 are smaller than the 33% statistical uncertainty in the determination of within the mask. Figure 14 shows the histogram of within the above mask, as well as in subsets of the data with various cuts in p. It shows that most sky pixels with reliable have low values, and that most of these pixels have large polarization fractions, above p = 5%.

As can be seen in Figs. 12 and 13, a similar structure for appears in all estimates in the selection mask, clearly showing that these structures are not caused by a single subsection of the data. We note that, outside the mask, shows structures similar to those observed in the mask. However, significant differences appear in some regions, in particular between odd and even Planck surveys. We attribute those to an imperfect bandpass mismatch correction or to the fact that no ADC correction has been applied here. We have also conducted tests in order to quantify the possible noise-induced bias on . Those are described in Appendix C. Figure 15 shows the map of when the resulting estimate of the bias has been subtracted. Comparison with Fig. 12 shows that the effect of bias essentially reduces low values, but does not explain the patterns observed in the map. We therefore conclude that the structures seen in the map of the polarization angle dispersion function are real, rather than being induced by noise and/or bias. In the rest of the analysis carried out here, we use the map of derived from the full survey and only consider pixels where the S/N on as derived from our Monte Carlo analysis is larger than 3. The resulting map is shown in Fig. 12.

The maps of reveal that regions with large variations of the polarization direction are organized in an intricate network of elongated and filamentary structures, some spanning several tens of degrees in the sky. The filamentary aspects of the maps generally persist in regions of low . These structures seem to have little to no correlation with structures in the total intensity map, except for the few degrees along the inner Galactic plane that systematically show low values of . Further analysis of the angular distribution function and the comparison with the polarization fraction are presented in Sect. 4.3.

Figure 16 shows the values of the observed averaged in bins of p as a function of the lag value. As expected, the angle dispersion function increases steadily with increasing lag. Lower values of systematically correspond to higher p values, as discussed in Sect. 4.3. Figure 17 shows details of for a few selected regions.

4. Discussion

In this section we analyse the observed variations of the polarization fraction and angle at 353 GHz and discuss the possible implications in terms of dust physics and magnetic field structure.

4.1. Maximum polarization fraction

The maximum dust polarization fraction (p_max) is a parameter of importance for the physics of dust and its alignment with respect to the magnetic field, because it directly provides an upper limit to the intrinsic polarization fraction, p₀, for the optimal orientation of the magnetic field, i.e., in the plane of the sky. It is also important for the CMB component separation in polarization, as it sets the maximum amplitude of dust polarization. The observed p values are, however, affected by averaging in the beam and along any given LOS. Variations of the B direction within the beam or along the LOS necessarily result in a decrease of the observed p. Similarly, dilution by an additional unpolarized source of emission, such as free-free or spinning dust emission, can only decrease p. Therefore, derived values of p_max can only be lower limits to the intrinsic polarization fraction p₀.

Here, we use the Planck maps at 353 GHz to evaluate p_max. Because p is a biased quantity and noise depends upon the data resolution, the observed maximum polarization fraction depends upon resolution. It is therefore crucial to take uncertainties into account. Figure 18 shows the sky fraction, f_sky(p>p_v), where the observed polarization fraction is higher than a given value p_v as a function of that p_v. The various curves are for data resolutions of 1°, 30′, and 15′. The coloured area shown correspond to f_sky(p ± 4σ_p>p_v) for the various resolutions.

At low f_sky values and high resolutions, high values of p are observed. Inspection of the maps indicates that these are point-like objects, either isolated pixels or actual point sources. Because we are interested in diffuse emission only, these isolated values are ignored in evaluating p_max. Table 2 lists the maximum and median values of p at different resolutions. It also shows the maximum value of p − 4σ_p observed at each resolution. We use the average of these values as a conservative estimate of p_max and find p_max> 19.8%. This indicates that, in the most favourable conditions for dust alignment, the intrinsic polarization fraction p₀ is larger than 19.8%.

The value for p_max derived here is significantly larger than generally found in polarization measurements in emission. Apart from the large 10–20% values observed with the Archeopsexperiment (Benoît et al. 2004), most previous observations have reported polarization fractions of less than 10%. Physical interpretation of this high value in terms of the physics of dust alignment and optical properties is beyond the scope of this paper. Note that a detailed comparison between the polarization fraction seen in emission in the Planck data and that seen in extinction in the visible is presented in Planck Collaboration Int. XXI (2015).

4.2. Polarization fraction vs. column density

We now analyse the variations of the polarization fraction p with dust column density. We use the Bayesian mean posterior estimate of p described in Sect. 2.3 and shown in Fig. 4, computed at 1° resolution. For the dust optical depth map, we use the map of τ₃₅₃ derived in Planck Collaboration XI (2014) computed at 1° resolution and we adopt their conversion factor from τ₃₅₃ to gas column density, derived from the correlation with HI data: (12)Figure 19 shows the distribution of the observed polarization fraction p as a function of N_H, as derived from dust optical depth, both for the sky fraction shown in Fig. 4 and for the same region but excluding the inner Galactic plane (i.e., excluding ℓ_II < 90deg or ℓ_II> 270deg, | b_II | < 2°). As can be seen in the figure, the plot shows both considerable scatter at a given N_H, and also systematic trends with N_H. The scatter is remarkable: for all column densities below N_H = 10²² cm^-2 (A_V ≃ 5 mag⁸), the polarization fraction spans all the values between the noise limit and maximum values up to 15–20%. At higher column densities, there is a sharp decrease of the values but the scatter, down to the noise limit is still present. The sensitivity of the Planck polarization measurements allows for the first time to detect a behaviour which is more complex than a power-law dependence of p with N_H. The scatter in p may be due to depolarization caused by fluctuations of the magnetic field orientation along the LOS or in the beam, and/or by intrinsic variations in p. Possible origins of this scatter are analysed in Planck Collaboration Int. XX (2015).

The largest polarization fractions are reached in the range of column densities 2 × 10²⁰ cm^-2<N_H < 10²¹ cm^-2 (0.1 <A_V < 0.5 mag). We observe an ensemble average polarization fraction, ⟨ p ⟩ of 7% at N_H = 10²¹ cm^-2 (A_V = 0.5 mag). The average values of p at lower column densities are not discussed in this paper, because a proper treatment would require a careful analysis of the residual bias in the method used to derive p. This will be the subject of an upcoming paper. At larger N_H (10²¹ cm^-2<N_H < 1.5 × 10²² cm^-2, 0.5 <A_V < 8 mag), the bulk of the p values are below p ≃ 10% and the maximum values show a steady decline. Over that range of column densities, the average polarization fraction ⟨ p ⟩ decreases down to ≃4%. We observe a sharp drop in ⟨ p ⟩ starting at about N_H ≃ 1.5 × 10²² cm^-2 (A_V ≃ 8 mag). Above N_H ≃ 4 × 10²² cm^-2 (A_V ≃ 20 mag) values of p are systematically below 4% with an average value of ⟨ p ⟩ ≃ 1–2%.

Toward nearby dense cores (n_H> 3 × 10⁴ cm^-3, size ~0.1 pc, N_H> 10²² cm^-2) the polarization fraction is observed to decrease systematically with N_H. This effect contributes to the sharp drop observed at 2 × 10²² cm^-2 (Fig. 19, bottom panel). Inspection of the Planck polarization map at 353 GHz shows many examples of such dips in p associated with nearby dense clouds. A systematic statistical study in the vicinity of Planck cold clumps will be presented in a forthcoming paper. Figure 20 shows the example of the dark cloud L134 (Tucker et al. 1976; Mattila et al. 1979) which is located at high Galactic latitude in the otherwise highly polarized Aquila Rift. L134 is one of the coldest Cold Clumps in the Planck catalogue (Planck Collaboration XXIII 2011). It is clearly seen that p can be as large as 10% in the external regions and decreases to values as low as 1% at the column density peak. This behaviour appears to be common in the high-latitude sky and confirms previous studies. Such a decrease of the polarization fraction toward large column densities on small (~a few 0.1 pc) scales was reported previously in ground-based measurements of polarization both in emission (Ward-Thompson et al. 2000; Matthews & Wilson 2000) and extinction (e.g., Gerakines & Whittet 1995; Whittet et al. 2008). This is usually interpreted as being due to a gradual loss of alignment of dust grains in dense shielded regions. In the likely hypothesis that dust alignment processes involve radiative torques responsible for the rotation of dust grains (Draine & Weingartner 1996; Hoang & Lazarian 2008), polarization in externally heated clouds is expected to drop off in the most shielded regions. The sharp decrease of p observed for N_H> 10²² cm^-2 in Fig. 19 is roughly consistent with such a scenario.

However, an increase of column density in the Galaxy is not necessarily associated with an increase of shielding, and the decrease of the polarization fraction with increasing column density could also be due to fluctuations in the orientation of the magnetic field along a long LOS, causing depolarization. In order to shed light on this depolarization effect, the companion paper Planck Collaboration Int. XX (2015) compares the polarization properties of the Planck dust emission with maps of polarized emission computed in simulations of MHD turbulence. The simulations are anisotropic to allow for an analysis of the influence of a large-scale magnetic field combined with a turbulent field. The polarized dust emission is computed using a uniform dust intrinsic polarization fraction p₀ = 20%. A large scatter in the polarization fraction p per bin of column density and a decrease of the maximum (and mean) values of p with N_H are found in the simulated maps, similar to those observed. The analysis reported in Planck Collaboration Int. XX (2015) however does not encompass the specific case of dense cores, but we cannot rule out that fluctuations of the magnetic field direction also contribute to depolarization within dense cores.

As shown in Fig. 19, which displays the dependence of p on N_H over the intermediate Galactic latitude sky and in regions excluding the inner Galactic plane, most lines of sight with very low p values are within the inner Galactic plane. The large gas column densities in the inner Galaxy (A_V> 20 mag) arise both in massive star forming regions (i.e., dense gas with n_H> 3 × 10⁴ cm^-3 for regions around 0.3 to 1 pc), but also along long lines of sight (say a few kpc) sampling mostly low density gas in the Molecular Ring. We argue that the contribution from such star forming regions in the inner Galaxy is small in the Planck maps at a resolution of 1°, because such regions have angular sizes smaller than 1′ if they are located further than 2 kpc from the Sun. The tail of high column densities in the inner Galaxy is therefore mostly due to long lines of sight sampling low density gas.

For lines of sight toward the inner Galactic plane, a related question is whether they are probing a dense cold medium, shielded from the ambient UV field, or if they result from the accumulation of low density material distributed over large distances. The apparent dust temperature can in principle be used to discriminate between these two situations. Figure 21 shows the distribution of the apparent dust temperature (T_obs), as derived from the dust SED fitting in Planck Collaboration XI (2014) using a modified grey-body fit, as a function of column density. As discussed in Planck Collaboration XI (2014) the apparent dust temperature steadily decreases with increasing column density, up to N_H ≃ 10²² cm^-2. The figure shows that, at higher column densities, T_obs increases again with N_H. The bulk of the large column densities above about 3 × 10²² cm^-2 therefore probes material in which dust is warmer than in the cold shielded cores, because it resides either in the low density medium, weakly shielded from the UV field of the inner Galaxy, or close to star-forming regions. In this case, the observed decrease of p is unlikely to be due to radiative transfer effects alone. This is taken as additional evidence that the fluctuations of the magnetic field direction could be the main origin of the decrease of p with column density.

4.3. Polarization fraction vs. angle dispersion function

Figure 12 shows the distribution of computed as described in Sect. 3.3 from the full survey at 353 GHz for 1° resolution and with a δ = 30′ lag used in the analysis.

The map of exhibits a wide range of values. A striking feature of the map is the existence of confined regions of high values, often reaching 50° to 70°, which are organized in an intricate network of filamentary structures, some of which span more than 30° in length. Figure 17 shows maps of selected regions around some of these high regions. Inspection of the polarization maps shows that these filamentary features generally lie at the boundary between regions with uniform, but different, magnetic field orientations on the sky. In order to quantify this statement, we decompose the map using the watershed morphological operator (e.g., Beucher & Meyer 1993). This allows us to segment the maps into a set of connected cells, equivalent to adjacent catchment basins, separated by the larger filamentary structures. The cells defined in this way are shown in Fig. 22 (top) for the first Galactic quadrant.

In the decomposition the cells are separated by one map pixel, an amount chosen to match the resolution of the map so that the values in cell areas do not include the high values present in the filamentary structures. Cells defined in this way have an average surface area of a few square degrees. For each cell i we compute the average polarization angle and its dispersion σ_ψ,i. We then consider each pair of adjacent cells, computing the difference between their average polarization angles and its associated uncertainty from . We also compute the average of over the boundary between the two cells in the pair. Figure 22 shows the fraction f of cell pairs with plotted as a function of . It can be seen that about 70% of cell pairs are separated in average angle by more than the uncertainty. The calculation was repeated for the four Galactic quadrants and led to the same conclusion.

Maps computed at larger lags look similar to those shown in Fig. 12, although with wider filamentary features, due to the larger scale of the analysis. Maps computed at smaller lags show filamentary features at the same locations as in Fig. 12, which indicates that the structures are in general unresolved. We also derived maps of at higher resolution. However, the noise and bias on increase quickly at higher resolution, which makes it impossible to follow the structure of the filamentary features down to the full Planck resolution of 5′ in most regions of the sky.

Comparison of in Fig. 12 and the observed polarization fraction map of Fig. 4 on large scales shows that, overall, the filamentary features of high correspond to low values of p. A similar trend was observed previously in the OMC-2/3 molecular clouds regions by Poidevin et al. (2010), using 14′′ resolution polarimetry data at 353 GHz. The Planck maps show that this is a general trend, as confirmed by the plot in Fig. 23, which shows that p and are anti-correlated. The best-fit correlation shown is given by (13)with α = −0.834 and β = −0.504, where p is unitless and is in degrees. Low p regions often correspond to regions where the observed polarization direction ψ changes. This result is in line with the findings of the previous section and further supports the view that variations in the magnetic field orientation play an important role in lowering the observed polarization fraction, as a result of integration along the LOS and/or within the beam.

The above results are compared with those inferred from MHD simulations in Planck Collaboration Int. XX (2015). The simulations clearly show an anti-correlation between and p, with a slope similar to that observed in the data. It is worth noting that in the noiseless simulations, the observed trend cannot be produced by the bias on resulting from higher uncertainties in polarization angles in regions of low signal and/or polarization fraction. It results from averaging effects of the polarization angle along the LOS. In brief, fluctuations of the magnetic field direction weaken the apparent polarization fraction, especially when the large-scale field tends to be aligned with the LOS.

The regions of large bear a morphological resemblance to features detected in maps of radio polarized emission, so-called “depolarization canals” (e.g., Haverkorn et al. 2000) and regions of high polarization gradient (Gaensler et al. 2011). The radio depolarization canals arise from Faraday rotation effects: they are thought to be due to either differential Faraday rotation (and hence depolarization) within synchrotron emission regions or discontinuities in foreground Faraday rotation screens (e.g., Fletcher & Shukurov 2007). The observed positions of depolarization canals vary with radio frequency and do not correspond to true physical structures. On the other hand, regions of high radio polarization gradient are somewhat similar to our filamentary features of large , insofar as each can be traced back to physical discontinuities. However, both the physical quantities that undergo a discontinuity to produce the phenomenon (free-electron density and LOS field component for the former versus sky-projected field orientation for the latter) and the places where the discontinuities occur (foreground Faraday rotation screens for the former versus dust-emitting regions for the latter) are unrelated. Therefore, one does not expect any one-to-one correspondence in morphology. Indeed, this is what we observe when comparing the distribution of the Planck filamentary features with maps of constructed using the 1.4 GHz synchrotron maps of Wolleben et al. (2006) and Testori et al. (2008).

Nevertheless, their morphological resemblance is rooted in a fundamental property of magnetized turbulence that goes beyond the detailed nature of the flows and of the specific ISM phase. Sub- or trans-sonic, super-Alfvénic turbulence in the warm ionized medium (WIM) is found to create filaments with high polarization gradients similar to those observed in the radio. They are shown to result from vorticity or shock compression (Gaensler et al. 2011). In the neutral ISM, both compressible and incompressible non-ideal magnetized turbulence generate filamentary structures with high at small scales similar to those observed (Planck Collaboration Int. XX 2015; Falgarone et al. 2015; Momferratos 2014). In the incompressible case, these filamentary structures are shown to follow current sheets of high intensity.

Finally, Fig. 24 shows the distribution of with column density.Several results are visible in these plots: (1) there is a large scatter of at all column densities, in particular below N_H ≃ 2 × 10²² cm^-2; (2) the largest values of are reached at low column densities, a fact that is well illustrated in the maps of Fig. 17; and (3) there is no more dependence of with N_H than found in that of p with N_H: the mean value of barely decreases with N_H over the narrow range 5 × 10²¹<N_H < 2 × 10²² cm^-2, then it increases sharply up to N_H ≃ 4 × 10²² cm^-2 and decreases again at larger column densities in the plot that includes sightlines across the Galactic plane. The small range of column densities over which increases with N_H is the same as that over which the polarization fraction (p) drops sharply (Fig. 19). This range corresponds to long lines of sight across the Galactic disk at Galactic latitudes >2deg. We take this result as an additional indication that, in that range of column densities, the drop of p is due to the fluctuations of the magnetic field orientation along the LOS and/or within the beam. This does not preclude other explanations, such as reduced dust alignment, for the low polarization fractions observed toward more opaque lines of sight. The dust polarization in dense Planck cold clumps, such as those described in Planck Collaboration XXII (2011), will be analysed in a forthcoming publication.

4.4. Dust vs. synchrotron polarization

In this section we compare the dust polarization as seen at 353 GHz with the synchrotron polarization that dominates at much lower frequencies. Our aim is to test how much the complementary observables trace the same magnetic fields and how their polarization properties are affected by the irregular component of the field. These comparisons tell us not only about the fields but also about the relative distributions of dust grains and relativistic electrons.

When Faraday rotation is negligible, the synchrotron and dust emission are both linearly polarized perpendicular to the local sky-projected magnetic field. The emissivities, however, have different dependencies on the magnetic field strength: the dust emission does not depend on the field strength, whereas the synchrotron emissivity is given by , where n_e is the density of relativistic electrons and γ is the power-law index of the relativistic-electron energy spectrum (typically γ ≃ 3, so that ). Synchrotron and dust polarization are also affected by different depolarization mechanisms. Some differences arise when the magnetic field in dust clouds differs from that in the diffuse synchrotron-emitting medium. Other differences arise because of the emissivity dependence on the field strength that weights the emission differently along the LOS. Any single direction may have a combination of these effects. We would therefore expect to see similar polarization structures where the particles sample the same average field, but not identical structures. A correlation analysis between dust and synchrotron polarization is also reported by Planck Collaboration Int. XXII (2015). Their cross-correlation between the Q and U maps at WMAP and LFI frequencies with the corresponding maps at 353 GHz shows that some of the polarized synchrotron emission at intermediate Galactic latitudes is correlated with dust polarization in both polarized intensity and angle. We might further expect to see statistical correlations even where the irregular component perturbs the large-scale magnetic field, but the degree of the correlation is complicated to predict.

The data sets are described in Sect. 2.5. Figure 25 compares both the polarization fraction and the polarization angle of the dust emission at 353 GHz (the Bayesian estimates where the S/N of p is greater than 3) with polarized synchrotron emission at 30 GHz. The comparison between the polarization angles is straightforward, because synchrotron is dominant and there is little Faraday rotation at that frequency. The comparison between the polarization fractions is more complex, however, because in the microwave and radio data there are additional total intensity components, such as free-free and anomalous microwave emission. To avoid contamination from anomalous microwave emission at 30 GHz, we begin with the 408 MHz map of Haslam et al. (1982) for total intensity. We correct for free-free emission as described in Sect. 2.5. This correction is approximate, but the synchrotron component dominates at low frequencies. Then we extrapolate the corrected 408 MHz synchrotron total intensity to 30 GHz in order to construct the polarization fraction, assuming a spectral index of − 2.88 (see, e.g., Jaffe et al. 2011). Note that a change in the constant value adopted for this index simply shifts the synchrotron polarization fraction systematically up or down and does not affect whether there would be an observed correlation. However, any spatial variations in the index that are not accounted for remain a limitation of this simple approach; they constitute an effective noise term that may blur the correlations we look for below but would not produce a spurious correlation.

The left two columns of Fig. 25 show the Galactic plane (| b_II | < 5°), while the right two show the results for the rest of the sky. The correlations are quantified by linear fits and Pearson correlation coefficients (r) listed in Table 3. In all but the fourth quadrant, there is weak but visible correlation (r> 0.3) in the polarization fraction in the plane, where the polarized intensity is strong. In the plane, we also see that the polarization angles remain near zero, i.e., the LOS integrated “apparent” magnetic field remains largely parallel to the plane. This confirms that at the largest scales probed through the full disk, the synchrotron and dust average over roughly the same structured magnetic fields. With a few notable exceptions, however, there is little correlation away from the plane, where isolated local structures and the irregular field component become more important.

A most interesting region for comparison is the second quadrant containing the Fan. In the Galactic plane this shows a relatively strong correlation (r = 0.47) in polarization fraction (as does the third quadrant to a lesser degree). Out of the plane, the correlation in p disappears. But we still see a correlation in the polarization angles off the plane, where they both remain concentrated around zero, indicating that the apparent magnetic field is parallel to the plane even at latitudes above 5°.

A second interesting region for the comparison is the first quadrant, where the sky is dominated by the radio loop I, i.e., the North Polar Spur (NPS). Here the high-latitude polarization angles show correlation where the two observables clearly trace the same magnetic fields.

We also compared the dust polarization angle dispersion with the polarized synchrotron emission at 1.4 GHz where it is subject to significant Faraday rotation effects. In Fig. 1 of Burigana et al. (2006) (based on data described in Sect. 2.5), the polarization fraction shows strong depolarization of the synchrotron emission within 30° of the plane, with the exception of the Fan region in the second quadrant. Much of the depolarization is so-called “beam depolarization”. A diffuse background source viewed through the roughly 1° beam results in emission co-added along slightly different lines of sight that pass through different turbulent cells; polarized emission even with initially uniform polarization angles gets Faraday-rotated differently and cancels out. One might then expect that the resulting synchrotron polarization fraction would anti-correlate with the dust polarization angle dispersion. Lines of sight toward highly turbulent regions should have low synchrotron polarization due to Faraday effects and high dust polarization angular dispersion. Such correlations are not generally apparent; however, in some regions such as the second quadrant dominated by the Fan we see this effect, implying that the dust and synchrotron emission in the Fan are tracing some of the same turbulent magnetic fields.

Finally, it is instructive to compare the dust polarization fraction with Faraday RMs of extragalactic radio sources. Recall that RMs are proportional to the LOS field component (which is positive/negative if the field points toward/away from the observer) times the free-electron density and integrated along the LOS, whereas the dust polarization fraction is an increasing function of the inclination angle of the magnetic field to the LOS. Therefore, if the large-scale field is coherent along the LOS through the Galaxy, then a field orientation close to the LOS tends to make RMs of extragalactic sources large (in absolute value) and the dust or synchrotron polarization fraction small, whereas a field orientation close to the plane of the sky does the opposite. As a result, one might expect a rough anti-correlation between RMs of extragalactic sources and dust polarization fraction. However, only a very loose anti-correlation may be expected at best because: (1) Faraday rotation and dust emission take place in different environments, with possibly different magnetic field directions; (2) RMs depend not only on the field inclination to the LOS, but also on the total field strength and on the free-electron column density; and (3) the LOS field component could undergo reversals, which would decrease RMs without correspondingly increasing the dust polarization fraction. Similarly, one might expect a rough positive correlation between RMs of extragalactic sources (again in absolute value) and dust polarization angle dispersion, because if the large-scale field is globally oriented closer to the LOS, the dust polarization angle is more sensitive to the fluctuating field.

Using the catalogues of Brown et al. (2003), Brown et al. (2007), Taylor et al. (2009), and Van Eck et al. (2011), in Fig. 26 we compare the RM of each source with the properties S and p of the polarized emission in the corresponding map pixel. This figure confirms the expected trends for the large-scale field in the Galactic plane, using both synchrotron and dust emission. Away from the plane (not shown) where more local structures dominate, we find no correlations. Because of the considerations outlined above, further work modeling the magnetized ISM on large scales and studying individual regions in detail will be necessary to understand the nature of the relationships among these observables.

We compared the dust observables with the standard deviation σ_RM of the extragalactic RM measures binned into low-resolution (N_side = 16 or ~ 4°) pixels. Again, in the plane, we found the expected loose anti-correlation (r = −0.3), where regions of high RM variation have low dust polarization fraction. In this case, the anti-correlation is likely due to changes in the field orientation, where lines of sight toward more turbulent regions measure both increased RM variation and decreased dust polarization due to depolarization effects. In the framework of this interpretation, we expect a weak positive correlation between the σ_RM and the dust angle dispersion function, an effect that we barely see as a hint of a correlation, with r = 0.2.

Lastly, we compared the degree of polarization of the extragalactic sources themselves with that of the diffuse Galactic dust emission in the same direction. Because a source’s degree of polarization anti-correlates with the standard deviation of the RM due to turbulent depolarization (Haverkorn et al. 2008), we might expect a positive correlation in the degree of polarization of the sources compared to the Galactic dust emission if the two observables are affected by the same turbulence. We find a hint of this correlation (r = 0.2) in the plane, but not at higher latitude.

The Planck polarization data at 353 GHz provide a new tracer of magnetic fields and an important complement to radio observations due to the different origins of the photons. This first look at the comparison of these observables confirms the expected large-scale correspondence as well as interesting correlations in the Fan and NPS regions. We find only weak correlations over much of the sky where the effects of local structures and the irregular field component dominate. This is not surprising but is nonetheless worth noting as it highlights the importance of, as well as the challenges inherent to, combining these data to build a coherent picture of the Galactic magnetic fields.

5. Conclusions

We have presented the Planck large-scale maps of polarized thermal dust emission at 353 GHz. This emission arises from non-spherical grains aligned with respect to the magnetic field. These data allow us for the first time to study dust polarization over large angular scales and open the field for many detailed studies to come.

The dust polarization fraction p displays a large scatter at all column densities below N_H ~ 10²² cm^-2. The maximum p is high, and we derive a lower limit to its maximal value of p_max = 19.8%. The highest polarization fractions are observed in a handful of individual regions, usually located in intermediate to low column density parts of the sky.

The large-scale spatial distribution of p shows a modulation of its upper values that is in general agreement with predictions of the general magnetic field structure of the MW, as constrained previously from synchrotron and RM data.

In addition to the large scatter of p, from the noise limit to more than 15%, there is a tendency for both its ensemble average and maximum value at a given N_H to steadily decrease with total column density. However, the decrease is shallow below N_H = 10²² cm^-2, becoming steeper between N_H = 10²² cm^-2 and N_H = 4 × 10²² cm^-2, and reaches a somewhat constant value of 2% above N_H = 4 × 10²² cm^-2.

The Planck polarization data at 353 GHz also allow precise measurements of the polarization direction ψ over most of the sky. Rotated by 90°, this direction shows the orientation of the apparent magnetic field projected on the plane of the sky, weighted by the dust emission and integrated on the LOS. The polarization angles in the Galactic plane are observed to be consistent with B lying mostly in the plane, as strongly suggested by previous synchrotron measurements. This is particularly true in the inner MW and in the highly polarized Fan region.

In order to characterize the structure of this apparent field, we compute a local measure of the dispersion of polarization angles at a given lag, the polarization angle dispersion function . It increases with lag, as previous observations have shown at smaller scales in specific regions. The sky distribution of reveals a spectacular network of unresolved filamentary structures with large values. This is the first time such structures have been observed for dust polarization. When they can be followed down to the Planck resolution, their widths are smaller than the beam, and some of them span large angular distances (several tens of degrees). These filamentary structures anti-correlate with p, in the sense that regions with maximal angle dispersions correspond to the lowest polarization fractions. We show that, in the large fraction of the sky we study, this is not due to any noise-induced bias on and is therefore a real effect. The filamentary structures appear to be separating regions that have different but homogeneous field orientations.

The anti-correlation between and p suggests that fluctuations of the magnetic field orientation have a major contribution to depolarization. This is also found in simulations of anisotropic MHD turbulence, without variations of the alignment properties of dust grains (see details in Planck Collaboration Int. XX 2015).

The filamentary structures of high bear some resemblance to the depolarization canals that are observed at radio frequencies and attributed to Faraday rotation effects, although there is no correspondence at small scales, which comes as no surprise because they have a different origin.

We compared the dust polarization fraction and angle with polarized synchrotron data. There are indications that the two tracers see the same magnetic field orientation, particularly interesting to see in the Fan region and the North Polar Spur, but that the detailed distributions of dust and high-energy electrons must be different in order to explain the observed maps. We infer a loose statistical correlation between extragalactic-source RMs and both the dust polarization fraction p and the angle dispersion function (negative and positive, respectively).

---

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, J.A., 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/M CTES (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. 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 and from a joint agreement between University of São Paulo, Brazil, and COFECUB, France (grant nos. USP 2007.1.433.14.2 and COFECUB Uc Te 114/08). We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. Some of the results in this paper have been derived using the HEALPix package.

References

Online material

Appendix A: Noise estimates for Planck smoothed maps

Here, we show how to smooth polarization maps and derive the covariance matrix associated to the smoothed maps.

Appendix A.1: Analytical expressions for smoothing maps of the Stokes parameters and noise covariance matrices

Smoothing total intensity maps is straightforward, but this is not the case for polarization maps. Because the polarization frame follows the sky coordinates and rotates from one centre pixel to a neighbouring pixel whose polarization will be included in the smoothing, in principle the (Q,U) doublet must be also rotated at the same time (e.g., Keegstra et al. 1997). The issue is similar for evaluating the effects of smoothing on the 3 × 3 noise covariance matrix, though with mathematically distinct results. In this Appendix, we present an exact analytical solution to the local smoothing of maps of the Stokes I, Q, and U, as well as the effects of smoothing on their corresponding noise covariance matrix.

Appendix A.1.1: Smoothing of Stokes parameters

Figure A.1 presents the geometry of the problem. Let us consider a HEALPix pixel j at point J on the celestial sphere, with spherical coordinates (ϕ_⋆,θ_⋆). To perform smoothing around this position with a Gaussian beam with standard deviation σ_{1 / 2} = FWHM/ 2.35 centred at the position of this pixel we select all HEALPix pixels that fall within 5 times the FWHM of the smoothing beam (this footprint is sufficient for all practical purposes). Let k be one such pixel, centred at the point K with coordinates (ϕ_k, θ_k), at angular distance β from J defined by (A.1)The normalized Gaussian weight is then (A.2)and ∑ _kw_k = 1. Before averaging in the Gaussian beam, we need to rotate the polarization reference frame in K so as align it with that in J. For that the reference frame is first rotated by ψ_k into the great circle running through K and J, then translated to J, and finally rotated through −ψ_⋆. The net rotation angle of the reference frame from point K to point J is then (A.3)Due to the cylindrical symmetry around axis z, evaluating does not depend on the longitudes ϕ_⋆ and ϕ_k taken separately, but only on their difference (A.4)Using spherical trigonometry in Fig. A.1 with the HEALPix convention for angles ψ_⋆ and ψ_k, we find: To derive ψ_k and ψ_⋆ we use the two-parameter arctan function that resolves the π ambiguity in angles: (A.9)Because of the tan implicitly used, sinβ (a positive quantity) is eliminated in the evaluation of ψ_⋆, ψ_k, and .

We can now proceed to the rotation. It is equivalent to rotate the polarization frame at point K by the angle , or to rotate the data triplet (I_k, Q_k, U_k) at point K by an angle around the axis I. The latter is done with the rotation matrix (e.g., Tegmark & de Oliveira-Costa 2001) (A.10)Finally, the smoothed I, Q, and U maps are calculated by: (A.11)

Appendix A.1.2: Computing the noise covariance matrix for smoothed polarization maps

We want to compute the noise covariance matrix [ C] _⋆ at the position of a HEALPix pixel j for the smoothed polarization maps, given the noise covariance matrix [ C ] at the higher resolution of the original data. We will assume that the noise in different pixels is uncorrelated. From the given covariance matrix [ C ] _k at any pixel k we can produce random realizations of Gaussian noise through the Cholesky decomposition of the covariance matrix: where in the decomposition is the transpose of the matrix [ L ] _k and (G)_k = (G_I,G_Q,G_U)_k is a vector of normal Gaussian variables for I, Q, and U.

Applying Eq. (A.11) to the Gaussian noise realization, we obtain (A.14)The covariance matrix of the smoothed at the position J is given by (A.15)If the noise in distinct pixels is independent, as assumed, then ⟨ (G)_k (G)_i ⟩ = δ_ki, the Kronecker symbol, and so (A.16)which can be computed easily with Eq. (A.10).

Developing each term of the matrix, we can see more explicitly how the smoothing mixes the different elements⁹ of the noise covariance matrix: where we note that and depend on j and k. The mixing of the different elements of the covariance matrix during the smoothing is due not to the smoothing itself, but to the rotation of the polarization frame within the smoothing beam.

Appendix A.1.3: Smoothing of the noise covariance matrix with a Monte Carlo approach

For the purpose of this paper, we obtained smoothed covariance matrices using a Monte Carlo approach.

We first generate correlated noise maps (n_l, n_Q, n_U) on I, Q, and U at the resolution of the data using (A.23)where G_l, G_Q, and G_U are Gaussian normalized random vectors and L is the Cholesky decomposition of the covariance matrix [ C ] defined in Eq. (A.12).

The above noise I, Q, an U maps are then smoothed to the requested resolution using the smoothing procedure of the HEALPix package. The noise maps are further resampled using the udgrade procedure of the HEALPix package, so that pixellization respects the Shannon theorem for the desired resolution. The smoothed covariance matrices for each sky pixel are then derived from the statistics of the smoothed noise maps. The Monte Carlo simulations have been performed using 1000 realizations.

Both the analytical and the Monte Carlo approaches have been estimated on the Planck data and shown to give equivalent results.

Appendix B: Debiasing methods

Because p is a quadratic function of the observed Stokes parameters (see Eq. (1)) it is affected by a positive bias in the presence of noise. The bias becomes dominant at low S/N. Below we describe briefly a few of the techniques that have been investigated in order to correct for this bias. For a full discussion of the various debiasing methods, see the introductions in Montier et al. (2015a,b) and references therein.

Appendix B.1: Conventional method (method 1)

This method is the conventional determination often used on extinction polarization data. It uses the internal variances provided with the Planck data, which includes the white noise estimate on the total intensity (C_II) as well as on the Q and U Stokes parameters (C_QQ and C_UU) and the off-diagonal terms of the noise covariance matrix (C_IQ,C_IU,C_QU).

The debiased p² values are computed using (B.1)where is the variance of p computed from the observed Stokes parameters and the associated variances as follows: (B.2)The uncertainty on ψ is given by (B.3)where σ_P is the uncertainty on the polarized intensity: (B.4)In the case where I is supposed to be perfectly known, C_II = C_IQ = C_IU = 0, (B.5)Because it is based on derivatives around the true value of the I, Q, and U parameters, this is only valid in the high S/N regime. The conventional values of uncertainties derived above are compared to the ones obtained using the Bayesian approach in Fig. B.1.

Appendix B.2: Time cross-product method (method 2)

This method consists in computing cross products between estimates of Q and U with independent noise properties. In the case of Planck HFI, each sky pixel has been observed at least four times and the four independent surveys can be used for this purpose. Another option is to use half-ring maps which have been produced from the first and second halves of each ring. These methods have the disadvantage of using only part of the data, but the advantage of efficiently debiasing the data if the noise is effectively independent, without assumptions about the Q and U uncertainties.

In that case, can be computed as (B.6)where the sum is carried out either over independent survey maps or half-ring maps.

The uncertainty of p² can in turn be evaluated from the dispersion between pairs (B.7)

Appendix B.3: Bayesian methods (method3)

We use a method based on the one proposed by Quinn (2012) and extended to the more general case of an arbitrary covariance matrix by Montier et al. (2015a). We use the Mean Posterior Bayesian (MB) estimator described in Montier et al. (2015b). Unlike the conventional method presented in Sect. B.1, this method is in principle accurate at any signal-to-noise ratio. Figure B.1 compares the Bayesian predictions for p and ψ and their uncertainties σ_p and σ_ψ with those obtained from the conventional method (Eqs. (1), (2), and (B.1)–(B.4)) as predicted from the Planck data at 1deg resolution. As can be seen in the figure, the bias on p is generally important even at this low resolution. The conventional uncertainties are accurate only at low uncertainties, as expected because Eqs. (B.4) and (B.5) are obtained from Taylor expansion around the true values of the parameters. The difference in the uncertainties is the greatest for σ_ψ as the true value can only reach 52deg for purely random orientations.

Appendix C: Tests on the bias of

We have performed tests in which we used the Planck noise covariance matrices in order to check that the structures we observe in the maps of the polarization angle dispersion function are not caused by systematic noise bias. One of the tests (called ) consisted of assigning each pixel a random polarization angle ψ. The second one (called ) consisted of setting ψ to a constant value over the whole sky map, which leads to (except near the poles). In that case, changing ψ in the data was done while preserving the value of p and σ_p computed as in Eq. (B.2), through the appropriate modification of I, Q, and U. The tests also use the noise covariance of the data, so that the tests are performed with the same sky distribution of the polarization S/N as in the data. This is critical for investigating the spatial distribution of the noise-induced bias on . In both tests, we added correlated noise on I, Q, and U using the actual noise covariance matrix at each pixel, and computed the map of using Eq. (6) and the same lag value as for the Planck data.

Figure C.1 (upper) shows the histograms of the values obtained for these two tests, both for the whole sky and in the mask used in the analysis of the real data. histograms peak at the value of for Gaussian noise only (no signal, ). The corresponding map of does not exhibit the filamentary structure of the actual data shown in Fig. 12. Similarly, the test histograms of do not resemble that of the real data shown in Fig. 14.

The test is important for assessing the amplitude of the noise-induced bias, as Monte Carlo simulations show that assuming a true value of maximizes the bias. We therefore use this test as a determination of the upper limit for the bias given the polarization fraction and noise properties of the data. Figure C.1 (upper) shows that the histograms of the recovered values peak at . The histogram is also narrower in the high S/N region than over larger sky regions at lower S/N. In the high S/N mask, 60% of the data points have a noise-induced bias smaller than 1.6°, and 97% have a bias smaller than 9.6°. The maps of the bias computed for this test show a correlation with the map of . However, as shown in Fig. C.1 (lower panel), the effect of the bias (the size of the offset) is small at high values of for most pixels and can reach up to 50% for a larger fraction of points at lower values (say below ). This

bias correction does not significantly change the structure of the map shown in Fig. 12 and so, in particular, bias does not cause the filamentary structures observed. We note, however, that the noise-induced bias can change the slope of the correlation between and p.

All Tables

All Figures