© Planck Collaboration 2020

1. Introduction

Interstellar dust grains are heated by absorption of the interstellar radiation field (ISRF), the ambient ultraviolet (UV), visible, and near-infrared radiation produced by the ensemble of stars in the Galaxy. The grains cool via thermal emission, which is in the far-infrared/submillimetre, as determined by the equilibrium temperature corresponding to a balance between absorbed and emitted power. Thermal emission from the larger grains that dominate the mass in the grain size distribution can be modelled as that of a modified blackbody (MBB) with emissivity ϵ_ν = κ_νB_ν(T_d), where the absorption coefficient κ_ν depends on the dust properties (Kruegel 2003). The equilibrium temperature is observed to be of order 20 K (Planck Collaboration XI 2014) for the ISRF found in the bulk of the interstellar medium (ISM).

Starlight polarization, discovered by Hall (1949) and Hiltner (1949), was quickly ascribed to differential extinction by aspherical dust grains with a preferential alignment related to the configuration of the interstellar magnetic field (Davis & Greenstein 1949, 1951). Over the years, a number of theories have been put forward to explain how this alignment occurs and is sustained, despite gas collisions (see the review by Andersson et al. 2015). The mechanism favoured currently involves radiative torques acting on grains subject to anisotropic illumination (RAT; see, e.g., Hoang & Lazarian 2016).

For thermal processes, Kirchhoff’s law states that differential extinction implies differential emission and so the submillimetre thermal emission from dust grains is also polarized, orthogonally to that of extinction. Thus, for dust grains aligned with respect to the Galactic magnetic field (GMF), the observed emission is also partially linearly polarized (Stein 1966; Hildebrand 1988). Because the spin axis of a dust particle is perpendicular to its long axis and alignment is statistically parallel to the local orientation of the magnetic field, the polarization of starlight transmitted through interstellar dust reveals the average orientation of the magnetic field projected on the plane of the sky, whereas the direction of polarized emission is rotated by 90° with respect to the magnetic field.

Observations of this submillimetre emission from Galactic dust, in both total intensity and polarization, have drawn strong attention, thanks to the Planck¹ full-sky maps, whose sensitivity and sky-coverage largely supersede the previously-available data from ground-based, balloon-borne (e.g., de Bernardis et al. 1999; Benoît et al. 2004), and space observations (e.g., Gold et al. 2011).

Over the course of four years (2009–2013), Planck surveyed the entire sky in nine frequency bands, from 30 GHz to 857 GHz, providing the best maps to date of the cosmic microwave emission, with unprecedented sensitivity, and angular resolutions varying from 30′ at 30 GHz to 4.​​′8 at 857 GHz. All but the two highest-frequency channels (545 GHz and 857 GHz) were sensitive to linear polarization of the observed radiation. In these seven bands, most of the polarized signal is of Galactic origin, with polarized synchrotron emission dominating at the low-frequency end of the spectrum, and polarized thermal emission from Galactic dust dominating at the high-frequency end. At 353 GHz, which is therefore the highest-frequency polarization-sensitive channel of Planck, polarized thermal dust emission is about two orders of magnitude stronger than the polarized cosmic microwave background (CMB; Planck Collaboration I 2016). It is therefore the channel we use to study this Galactic emission, and several Planck papers have already provided analyses of earlier releases of this data to investigate the link between dust polarization and physical properties of the ISM, most notably the structure of the Galactic magnetic field, properties of dust grains, and interstellar turbulence. In Appendix A, we provide a summary of the main results of these Planck papers, to serve as a useful reference.

In this paper, one in a series associated with the 2018 release of data from the Planck mission (Planck Collaboration I 2020), we use all-sky maps of dust polarized emission produced from this third public release of Planck data (hereafter the Planck 2018 data release or PR3) to expand on some of these studies of the polarized thermal emission from Galactic dust. More specifically, our analysis first focuses on a refined statistical analysis of the dust emission’s polarization fraction and polarization angle over the full sky, in the fashion of Planck Collaboration Int. XIX (2015) but based on a post-processing of the Planck 2018 data that minimizes the contamination from components other than dust. One of the results from that paper, confirmed by a comparison with numerical simulations of magnetohydrodynamical (MHD) interstellar turbulence (Planck Collaboration Int. XX 2015), is the nearly inverse proportionality of the polarization fraction p and the local dispersion of polarization angles 𝒮. Here we propose an interpretation of this relationship, showing that it is a generic result of the turbulent nature of interstellar magnetic fields. We therefore further analyse the Planck data by considering the product 𝒮 × p, which allows us to search for deviations from this first-order relationship. Deviations might be related to changes in the properties of the dust or of its alignment with respect to the magnetic field. In the final part of the paper, we present an updated comparison of the dust polarized emission with stellar polarization data in the visible, following Planck Collaboration Int. XXI (2015), but with a much larger sample of stellar polarization data. For aspects of the analysis of polarized thermal dust emission related to component-separation, i.e., the angular power spectra and spectral energy distributions (SEDs) of the E and B modes, we refer the reader to Planck Collaboration XI (2020).

The paper is organized as follows. In Sect. 2, we present the Planck maps of Stokes parameters that are used in the subsequent analysis. In Sect. 3, we present the full-sky maps of thermal dust polarization derived from these Stokes maps. In Sect. 4, we present a statistical overview of these dust polarization maps over the full sky, using the tools and analysis introduced in Planck Collaboration Int. XIX (2015). In Sect. 5 we expand on this statistical analysis, looking for trends beyond the first-order correlations exhibited by the data. In Sect. 6, we update our comparison with the stellar polarization data, greatly expanding on the sample presented initially in Planck Collaboration Int. XXI (2015). Finally, Sect. 7 presents our conclusions. Seven appendices complete the paper. In Appendix A, as already mentioned, we offer a summary of the main results of earlier Planck papers dealing with the polarized thermal emission from Galactic dust. In Appendix B, we show complementary, variable resolution Stokes maps at 353 GHz and present the Stokes covariance maps that are used to assess the statistical uncertainties affecting Planck polarization data presented in this work. Appendix C describes our approach to estimating the systematic uncertainties in the data, based on a set of end-to-end (E2E) simulations. Appendix D explains the relationship of the polarization angle dispersion function 𝒮 to the polarization gradients commonly used in polarization studies at lower frequencies. Appendix E provides supplementary figures showing how the behaviour of polarization fraction with total gas column density is affected by the uncertainty on the Galactic zero level. Appendix F provides a demonstration of the inverse relationship between the polarization fraction p and the polarization angle dispersion function 𝒮, based on a phenomenological model of magnetized interstellar turbulence. Finally, Appendix G assesses the noise and systematics that affect the data used in the comparison of visible and submillimetre polarization properties (Sect. 6).

2. Processing Planck maps for Galactic science

The Stokes I, Q, and U maps at 353 GHz that we use in this paper are based on products from the Planck 2018 data release. The processing steps applied to the data to compute the Planck 2018 frequency maps are presented in Planck Collaboration II (2020) and Planck Collaboration III (2020) for the Low Frequency Instrument (LFI) and High Frequency Instrument (HFI), respectively. For HFI, the Q and U products used at 353 GHz make use of the polarization-sensitive bolometers (PSBs) only, ignoring the spider-web bolometer (SWB) data, as recommended in Planck Collaboration III (2020), while the rest, including I at 353 GHz, make use of the complete data set (PSB+SWB).

For our Galactic science applications, we use maps that result from post-processing with the Generalized Needlet Internal Linear Combination (GNILC) algorithm, developed by Remazeilles et al. (2011); this filters out the cosmic infrared background (CIB) anisotropies, a key feature for Galactic science. These GNILC maps, derived from the Planck 2018 maps, are presented and characterized in Planck Collaboration IV (2020), and so we simply recall a few key properties of this post-processing step in the next subsection (Sect. 2.1). The GNILC maps used here have a uniform resolution of 80′.

In Sects. 5.3, 5.4, and 6, where we require data at a uniform resolution that is finer than 80′, and where we are less concerned by the presence of CIB anisotropies, we use maps derived more directly from the Planck 2018 353 GHz Stokes maps and their covariance maps. The required post-processing to produce these alternative Stokes maps (ASMs) is also described below.

As another important post-processing step, we need to establish the desired zero level in the total intensity maps for Galactic dust emission, as described in Sect. 2.2.

2.1. GNILC and ASM post-processing

GNILC is a wavelet-based component-separation method that makes use of both spectral and spatial information to disentangle multidimensional components of the sky emission. In practice it combines data from the different Planck bands and outputs maps at any desired frequency.

In Planck Collaboration Int. XLVIII (2016), GNILC was applied to Planck 2015 total intensity data, effectively separating Galactic thermal dust emission and CIB anisotropies over the entire sky, while simultaneously filtering out noise and CMB contributions. In regions of low dust column density, it was found that the CIB anisotropies are well above the noise, correlated spatially, and provide a significant contribution to the emission power spectrum. We are interested in polarization properties for Galactic dust emission over the full sky, including high-latitude diffuse lines of sight, for which GNILC-processing significantly reduces contamination of the I map by CIB anisotropies.

For the Planck 2018 data release we go further, applying GNILC not only in total intensity, but also in polarization, thus providing maps of polarized Galactic thermal dust emission in which the contamination by polarized CMB emission and instrumental noise has been reduced.

The GNILC algorithm optimizes the component separation given the local variations of the contamination. At high Galactic latitudes and small angular scales, the local dimension of the Galactic signal subspace estimated by GNILC, i.e., the number of components in the Galactic signal, can be null because in this regime the data become compatible with a mixture of CIB, CMB, and noise². Therefore, the effective resolution of the GNILC dust maps is not uniform but variable over the sky, with an effective beam whose full-width at half-maximum (FWHM) increases from the Galactic plane towards high latitudes. The local resolution depends on the local signal-to-nuisance ratio, which varies differently for intensity and for E- and B-mode polarization³. Therefore, the optimal GNILC resolution should, by design, be different for total intensity and for polarization maps. However, for consistency in the astrophysical study of dust intensity and polarization, where the polarization fraction p = P/I is of interest, we adopt a common resolution by imposing that the variable resolution of the GNILC dust maps should be driven by the more stringent signal-to-nuisance ratio of the B-mode data. In practice, in the Galactic plane the signal-to-noise ratio (S/N) in polarization is sufficiently large to allow for the use of the nominal Planck resolution at 353 GHz, while for high Galactic latitude regions data smoothed to 80′ are required.

The GNILC method is also able to provide Stokes maps at a uniform resolution of 80′ over the entire sky, enabling the analysis of polarization properties over the entire sky at a common resolution. It should be noted that in this case, and to avoid oversampling, the output maps are subsequently downgraded from the original HEALPix⁴ (Górski et al. 2005) resolution N_side = 2048 to N_side = 128.

The equivalent ASM post-processing step is to subtract the total intensity CMB SMICA map (Cardoso et al. 2008; Planck Collaboration IV 2020) from the Planck 2018 total intensity map at 353 GHz. No subtraction of CIB anisotropies is performed. Compared to the dust signal at 353 GHz, the CMB polarized signal is small, less than 1% (Planck Collaboration IV 2020), and subtracting that would add noise unnecessarily.

2.2. Zero level for total intensity of Galactic thermal dust emission

We recall that Planck had very little sensitivity to the absolute level of emission and so the zero level of the maps of I must be set using ancillary data. This is of central interest for our study, because for the most diffuse lines of sight it directly impacts polarization fractions through p = P/I.

Planck 2018 HFI frequency maps, as delivered (Planck Collaboration III 2020), deliberately include a model of the CIB monopole. As a first step towards maps suitable for Galactic science, this needs to be subtracted. GNILC post-processing does not adjust the monopoles contained in the input maps and so the CIB monopole needs to be subtracted explicitly, frequency by frequency, as for ASMs. At 353 GHz the intensity of the model CIB monopole is 0.13 MJy sr⁻¹, or 452 μK_CMB using the unit conversion $287.5\mathrm{MJy}{\mathrm{sr}}^{-1}{\mathrm{K}}_{\mathrm{CMB}}^{-1}$ given in Planck Collaboration III (2020).

This CIB-subtracted total intensity map has a zero level that by construction is based on a correlation of the emission at high Galactic latitudes with the column density of the ISM traced by the 21-cm emission of H I at low column densities. Nevertheless, this Galactic offset needs to be refined. A favoured method is again based on a correlation of dust emission with H I, as described in Planck Collaboration VIII (2014), Planck Collaboration XI (2014), Planck Collaboration Int. XLVIII (2016), and Planck Collaboration III (2020). After the GNILC processing, we apply the same H I correlation procedure to the output maps of I, in particular finding that a Galactic H I offset of 36 μK_CMB should be added to the 353-GHz GNILC total intensity map used for polarization at the uniform 80′ resolution. The statistical error of about 2 μK_CMB is small compared to the systematic uncertainties that we now discuss.

Because the dust total intensity versus H I correlation has an upward curvature, the estimates of the offset and slope are dependent on the column density range used for the fit. Furthermore, there is an additional source of uncertainty, related to the possibly significant emission from dust that is in the warm ionized medium (WIM), and therefore associated with H II rather than with neutral hydrogen H I. The fractional contribution might be most important at low H I column densities, i.e., in the diffuse ISM.

To assess the systematic effect related to the WIM-associated dust, we rely on an estimate of the total column density of the WIM towards high Galactic latitudes by Gaensler et al. (2008), N_H, WIM = 8 × 10¹⁹ cm⁻². Assuming the same SED in the submillimetre per proton as per H atom, and using the results of Planck Collaboration Int. XVII (2014), this translates to 54 μK_CMB at 353 GHz. If all of the dust emission associated with the WIM were uncorrelated with the H I-associated dust, then this value would need to be added to the Galactic H I offset. On the other hand, part of any dust emission associated with the WIM is probably correlated with H I as well, and in the extreme case of 100% correlation, there would be no correction due to the WIM.

To account for this effect, we adopt a central value of 27 μK_CMB which, when added to the Galactic H I offset, gives a fiducial total Galactic offset of 63 μK_CMB (corresponding to 0.0181 MJy sr⁻¹ at 353 GHz), to be added back to the GNILC total intensity map at 353 GHz, after the CIB monopole subtraction. This fiducial value will be used in the rest of our analysis. It has an uncertainty that we estimate to be 40 μK_CMB (corresponding to ±0.0115 MJy sr⁻¹). As mentioned above, the offset affects the statistics of the polarization fraction of dust polarized emission. To quantify the effect of an offset uncertainty in the range estimated, we also use intensity maps resulting from the addition of a total Galactic offset of 23 μK_CMB (0.0066 MJy sr⁻¹) and 103 μK_CMB (0.0296 MJy sr⁻¹), referred to as low and high, respectively. Note, however, that these correspond to fainter and brighter intensity maps, leading to higher and lower polarization fractions, respectively.

The procedure to adjust the ASM intensity map at 353 GHz after CIB-monopole subtraction is the same. In this case the fiducial Galactic offset is 68 μK_CMB, a value that is, not surprisingly, very close to that for GNILC.

2.3. GNILC Stokes maps

For the 353-GHz data used here, after the adjustments of the zero level of I just discussed, the GNILC Stokes I, Q, and U maps are converted to astrophysical units using the already mentioned conversion factor $287.5\mathrm{MJy}{\mathrm{sr}}^{-1}{\mathrm{K}}_{\mathrm{CMB}}^{-1}$. The resulting GNILC Stokes maps at 353 GHz and uniform 80′ resolution are shown⁵ in Fig. 1. The total intensity map corresponds to the fiducial offset value. The GNILC Stokes maps at 353 GHz and variable resolution over the sky are shown in Appendix B, alongside the GNILC-processed covariance maps σ_II, σ_IQ, σ_IU, σ_QQ, σ_QU, and σ_UU that are used in Sect. 3.2 to estimate the statistical uncertainties on the dust polarization properties.

We note that for studies involving the polarization angle dispersion function 𝒮 (Sect. 3.3), we use Stokes maps and covariance maps that are further smoothed to a 160′ FWHM uniform resolution, and downgraded to N_side = 64.

2.4. Alternative Stokes maps (ASMs)

For ASMs, as a final step after converting to astrophysical units, we smooth the Stokes I, Q, and U maps uniformly to 10′, 20′, 40′, 60′, 80′, and 160′, downgrading the HEALPix resolution to N_side = 1024, N_side = 512, N_side = 256, N_side = 128, and N_side = 64, respectively. The covariance matrix maps σ_II, σ_IQ, σ_IU, σ_QQ, σ_QU, and σ_UU are consistently smoothed from Planck 2018 data to the same resolutions using the procedure described in Appendix A of Planck Collaboration Int. XIX (2015).

3. Full-sky thermal dust polarization maps

In this section, we present the maps of Galactic thermal dust polarization over the full sky, derived from the GNILC-processed Stokes I, Q, and U maps at uniform 80′ resolution.

3.1. Polarization fraction and angle maps

From the GNILC maps of Stokes parameters I, Q, and U at 353 GHz, we build maps of the polarized intensity P, polarization fraction p, and polarization angle ψ. The convention used for the Stokes parameters in the Planck 2018 data release is to measure polarization angles from the direction of the Galactic north and positively towards Galactic west in accordance with the HEALPix convention used in cosmology (see Planck Collaboration 2018, for further discussion). However, as in Planck Collaboration Int. XIX (2015), we conform here to the IAU convention, polarization angles ψ being counted positively towards Galactic east, and so they are computed simply by changing the sign of Stokes U in the Planck 2018 data. Thus

$$\begin{array}{c}\hfill P=\sqrt{Q^2+U^2}p=\frac{P}{I}\psi =\frac{1}{2}\mathrm{atan}2(-U,Q),\end{array}$$(1)

where the two-argument function atan2(−U, Q) is used in place of atan(−U/Q) to avoid the π-ambiguity. Conversely, the Stokes parameters can be recovered from the total intensity, the polarization fraction, and the polarization angle via

$$\begin{array}{c}\hfill Q=pIcos\left(2\psi \right)U=-pIsin\left(2\psi \right).\end{array}$$(2)

The presence of noise in the Stokes maps can bias the estimates of P, p, and ψ (Montier et al. 2015a,b), so that naive estimators $\widehat{P}$, $\widehat{p}$, and $\widehat{\psi }$ computed using Eq. (1) directly on the noisy data do not adequately represent the true values at low S/N. Alternative estimators have been developed, most notably for the polarized intensity and the polarization fraction (the bias on the polarization angle is usually negligible). For the polarization fraction, we use the modified asymptotic (MAS) estimator introduced by Plaszczynski et al. (2014) and defined through

$$\begin{array}{c}\hfill p_{\mathrm{MAS}}=\widehat{p}-{\varsigma }^2\frac{1-e^{-{\widehat{p}}^2/{\varsigma }^2}}{2\widehat{p}},\end{array}$$(3)

where ς is a noise-bias parameter that depends on the geometrical properties of the (assumed Gaussian) 2-dimensional distribution of the noise in (Q, U) space, assuming a noise-free total intensity I. From the 353-GHz GNILC covariance matrices at the uniform 80′ resolution, we can compute this noise-bias parameter and find that ς² <  10⁻³ over the full sky, which shows that the debiasing performed by the MAS estimator is small.

Because the noise in total intensity is small, this is a reasonable approach that can also be used to provide a MAS estimate of the polarized intensity, P_MAS. For notational simplicity, we hereafter drop the subscript “MAS” and write p to mean p_MAS and P to mean P_MAS. For the GNILC-processed 353-GHz data at the uniform 80′ resolution, the resulting polarized intensity P map is shown in Fig. 1 (bottom row), while the polarization fraction p and the polarization angle ψ maps are shown in Fig. 2. We note that the total intensity offset used for these maps is the fiducial one. The choice of offset has an impact on p (as we will discuss in Sect. 4.1) but not on ψ or P.

The overall structure of the polarization fraction and angle maps is consistent with that found over a smaller fraction of the sky in Planck Collaboration Int. XIX (2015). We note in particular that the Galactic plane exhibits low polarization fractions, except towards the “Fan” region, near the anticentre, and that the structures seen in p do not generally correspond to structures in total intensity. The polarization angle map ψ shows that the magnetic field is essentially parallel to the Galactic plane at low Galactic latitudes |b|, and the large-scale patterns at higher latitudes can be broadly interpreted as arising from the projection of the local magnetic field in the Solar neighbourhood (Planck Collaboration Int. XLIV 2016; Alves et al. 2018).

3.2. Estimation of uncertainties

There are several types of uncertainties that need to be taken into account in our analysis of the dust polarization maps.

First, there is statistical noise, whose contribution to the uncertainties can be estimated using the covariance maps of the GNILC-processed data. This was evaluated by performing a set of Monte Carlo simulations of Stokes I, Q, and U maps, taking the GNILC maps as means of a multivariate normal distribution with covariances given by the GNILC covariance maps⁶. A set of 1000 simulations was computed; results for a set half this size do not change significantly, confirming that 1000 is sufficient. From these simulations we computed 1000 maps of not only p and ψ, but also other derived quantities, such as the polarization angle dispersion function (Sect. 3.3). As discussed in Sect. 4, these are instrumental in detecting any remaining bias (after using the MAS estimator), by investigating whether statistical properties (e.g., the histogram of p) computed on the GNILC maps shown in Fig. 2, are compatible with the ensemble average of the same properties computed on the Monte Carlo simulations. When, as expected, the quantities in the polarization maps are unbiased, the standard deviations of these maps, and of any derived quantity that we compute using the simulations, yield reliable statistical uncertainties.

Using this approach, we computed polarization fraction and polarization angle uncertainty maps σ_p and σ_ψ (shown in Fig. 2). These are actually very close to the ones obtained using Eqs. (B.2) and (B.3) of Planck Collaboration Int. XIX (2015), which are valid at sufficiently high S/N in polarization p/σ_p⁷. Figure 3 shows the polarization S/N map for the GNILC-processed data at 353 GHz and uniform 80′ resolution. At this resolution, p/σ_p >  3 over most of the sky and thus the estimate of the S/N is robust (Montier et al. 2015a)⁸.

The statistical absolute uncertainty on polarization fractions is at most 3%, and the statistical uncertainty on polarization angles is completely negligible, at less than 0.1°. Furthermore, based on the results of Montier et al. (2015a), we are confident that the polarization angle bias is less than 10% of this value. Indeed, at 80′ resolution, 99.9% of the sky pixels have an effective ellipticity below 1.25. This quantity characterizes the asymmetry between the noise distributions on Q and U maps in a rotated reference frame that cancels correlated noise between the two. Montier et al. (2015a) show that in this case the bias on the polarization angle is at most of order 7–8% of the statistical uncertainty σ_ψ.

Second, we need to estimate the impact of residual systematics arising from the Planck data processing. This is accomplished via a set of 100 end-to-end (E2E) simulations that take a model sky as input, simulate the timelines of the instrument taking into account all known systematics, and process these simulated timelines with the mapmaking pipeline described in Planck Collaboration III (2020). These E2E simulations are described in detail in Appendix A of Planck Collaboration XI (2020). The statistical comparison between the input and output polarization maps, which we discuss in Appendix C, shows that the absolute uncertainties from residual systematics are estimated to be ±0.5% on p and ±8° on ψ.

We note that these E2E simulations include realizations of random data noise and so already include part of the statistical uncertainty that is addressed by the Monte Carlo simulations based on the covariance matrices.

Finally, as already mentioned in Sect. 2.2, the quantitative analysis of p towards diffuse lines of sight depends strongly on the value of the Galactic offset used to set the zero level of total intensity for Galactic dust emission. To take this source of uncertainty into account, following the discussion in Sect. 2.2 we consider a fiducial case in which the Galactic offset is 63 μK_CMB and also consider a range of ±40 μK_CMB about this central value.

3.3. Polarization angle dispersion function

The polarization angle dispersion function 𝒮, introduced in Planck Collaboration Int. XIX (2015) quantifies the local (non-)uniformity of the polarization angle patterns on the sky by means of the local variance of the polarization angle map at a certain scale parameterized by a lag δ. It is defined as

$$\begin{array}{c}\hfill \mathcal{S}(\mathit{r},\delta )={\displaystyle \sqrt{\frac{1}{N}\sum _{i=1}^N{[\psi (\mathit{r}+{\mathit{\delta }}_i)-\psi (\mathit{r})]}^2},}\end{array}$$(4)

where the sum extends over the N pixels, indexed by i and located at positions r + δ_i, within an annulus centred on r and having inner and outer radii δ/2 and 3δ/2, respectively. Regions where the polarization angle tends to be uniform exhibit low values of 𝒮, while regions where the polarization patterns are more chaotic exhibit larger values, with $\mathcal{S}=\pi /\sqrt{12}\approx 52°$ when the polarization angles are completely uncorrelated spatially.

A map of 𝒮 at 60′ resolution and using a lag of 30′, based on Planck 2013 data, was shown over a restricted region of the sky in Planck Collaboration Int. XIX (2015). We can now present the 𝒮 map over the full sky, based on the GNILC-processed Planck 2018 data release at 353 GHz. Because 𝒮 is built from the polarization angle ψ, it is independent of the value chosen for the total intensity offset. However, when computed at uniform 80′ resolution and using a lag δ = 40′, 𝒮 is still significantly biased (see Sect. 4.1). For this reason, we use maps smoothed to 160′ and adopt a correspondingly larger lag δ = 80′⁹. This map is shown in the top panel of Fig. 4. We computed the statistical uncertainty ${\sigma }_{\mathcal{S}}^{\mathrm{MC}}$ using the Monte Carlo approach discussed in Sect. 3.2, but based on the Stokes maps smoothed to 160′ resolution. The map of ${\sigma }_{\mathcal{S}}^{\mathrm{MC}}$ is shown in the bottom panel of Fig. 4. Quite large values, up to 14°, are reached in some regions, but we will see in Sect. 4.1 that this is compatible with the noise in the data (see also Sect. 3.5).

3.4. Relationship of 𝒮 to alternative estimators

Synchrotron studies in the radio domain frequently use another estimator of the uniformity of polarization patterns, the polarization gradient introduced by Gaensler et al. (2011) and defined as

$$\begin{array}{c}\hfill |\mathrm{\nabla }P|=\sqrt{{\left(\frac{\partial Q}{\partial y}\right)}^2+{\left(\frac{\partial Q}{\partial z}\right)}^2+{\left(\frac{\partial U}{\partial y}\right)}^2+{\left(\frac{\partial U}{\partial z}\right)}^2},\end{array}$$(5)

where y and z refer to an orthogonal coordinate system on the plane of the sky. We show in Appendix D that, as far as the Planck thermal dust polarization data are concerned, |∇P| is strongly correlated with 𝒮, though not perfectly because of the contribution from the polarized intensity in |∇P|. This can be mitigated by considering an angular version of the polarization gradient defined as (Burkhart et al. 2012)

$$\begin{array}{c}\hfill |\mathrm{\nabla }\psi |=\sqrt{{\left[\frac{\partial (Q/P)}{\partial y}\right]}^2+{\left[\frac{\partial (Q/P)}{\partial z}\right]}^2+{\left[\frac{\partial (U/P)}{\partial y}\right]}^2+{\left[\frac{\partial (U/P)}{\partial z}\right]}^2},\end{array}$$(6)

which encodes only the angular content of the polarization¹⁰. In Appendix D, we show not only that |∇ψ| is better correlated with 𝒮 than |∇P| is, but also that this can be demonstrated analytically, with

$$\begin{array}{c}\hfill \mathcal{S}(\mathit{r},\delta )\approx \frac{\delta }{2\sqrt{2}}|\mathrm{\nabla }\psi |,\end{array}$$(7)

the linear dependence of 𝒮 on the lag being revealed simply through a first-order Taylor expansion. We do not use this estimator |∇ψ| in the rest of this paper, but note that in practice it might be easier to compute than 𝒮.

3.5. Noise and bias in 𝒮

An estimate of the variance of 𝒮 due to noise is (Planck Collaboration Int. XIX 2015; Alina et al. 2016):

$$\begin{array}{cc}\hfill {\sigma }_{\mathcal{S}}^2(\mathit{r},\delta )=& \frac{{\sigma }_{\psi }^2(\mathit{r})}{N^2{\mathcal{S}}^2}{{\displaystyle (\sum _{i=1}^N\psi (\mathit{r}+{\mathit{\delta }}_i)-\psi (\mathit{r}))2}}^{}\hfill \\ \hfill & +\frac{1}{N^2{\mathcal{S}}^2}{\displaystyle \sum _{i=1}^N{\sigma }_{\psi }^2(\mathit{r}+{\mathit{\delta }}_i){(\psi (\mathit{r}+{\mathit{\delta }}_i)-\psi (\mathit{r}))}^2.}\hfill \end{array}$$(8)

Just like for p, noise on Stokes parameters Q and U induces a bias on 𝒮. Unlike for p, however, this bias can be positive or negative, depending on whether the true value is, respectively, smaller or larger than the value $\pi /\sqrt{12}\approx 52{}^{°}$ obtained for fully random polarization angles (Alina et al. 2016). As prescribed by Hildebrand et al. (2009) and Planck Collaboration Int. XIX (2015), we use the following debiasing scheme

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{db}}=\{\begin{array}{cc}\sqrt{{\mathcal{S}}^2-{\sigma }_{\mathcal{S}}^2}\hfill & \text{if}\mathcal{S}>{\sigma }_{\mathcal{S}},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$(9)

This expression works well for S/N on 𝒮 larger than 3, which we ensure by smoothing the Stokes maps. For notational simplicity, in the rest of this paper, we write 𝒮 to mean the debiased value 𝒮_db of the polarization angle dispersion function.

4. Statistics of thermal dust polarization maps

In this section, we provide a statistical analysis of the quantities represented in the Galactic thermal dust polarization maps derived above. We start by discussing the distribution functions of p, ψ, and 𝒮. We then examine the joint distributions of p and total gas column density on the one hand, and of 𝒮 and p on the other hand. Finally, we look into how one striking feature of these maps, i.e., the inverse relationship between 𝒮 and p, is well reproduced by relatively simple Gaussian models of the Planck polarized sky.

4.1. Distribution functions for p, ψ, and 𝒮

4.1.1. Polarization fraction

The distribution function (DF, or histogram) for p over the full sky is shown in Fig. 5 The solid red curve is the histogram for the GNILC map of p for the fiducial offset in I, while the solid blue and green curves are the corresponding histograms for the low and high offsets, respectively. These clearly show the significant effect induced by the uncertainty on the total intensity offset. We note, however, that the polarization fractions observed reach at least 20% for any choice of the total intensity offset, setting strong constraints for dust models.

For comparison, the corresponding dashed coloured curves are the means of the DFs from the 1000 Monte Carlo simulations. Compared to the solid curves, there is only a small bias, shifting the DF towards higher p in the tail of the distribution. This is less pronounced for the green curves (high total intensity offset) because for this case the statistical changes in I are less important¹¹.

The coloured regions encompassing the mean histograms show the minimum and maximum values of the histogram for any given bin of p over the 1000 samples, i.e., the envelope within which all 1000 histograms lie. Lines defining the edges of the envelope would themselves not be distribution functions; however, they give an idea of the possible spread of the p histograms with varying noise realizations.

It is of interest, for dust models in particular, to estimate p_max, the largest value of p over the full sky. To estimate this and provide further quantification, we compute, for each of the total intensity offset values, the 90th, 95th, 98th, 99th, and 99.9th percentiles for the GNILC histogram from the data, which we write as h(p), and for each of the 1000 Monte Carlo histograms, which we write as h(p_i), with 1 ≤ i ≤ 1000. From the latter we calculate the mean and the standard deviation, which gives an estimate of the statistical uncertainty of p_max in a single realization, such as the data.

These numbers are given in Table 1, alongside the corresponding values for the average p map over the 1000 Monte Carlo realizations, which we write $h(\overline{p_i})$, and those for the mean histogram over the 1000 realizations¹², which we write $\overline{h(p_i)}$. The percentiles for the average p map are always very close to those for the data, which is not surprising because the data were taken as the mean for the Monte Carlo realizations. More interestingly, the percentiles h(p_i) are systematically larger than the corresponding values for the data, with very low statistical uncertainties. We note that this discrepancy is significantly smaller for the high total intensity offset than for the low total intensity offset, at least for the highest percentiles. This shows that p_max from the data is likely biased by a similar amount and is to be adjusted accordingly. We also point out that the percentiles for the mean histogram $\overline{h(p_i)}$ are larger still, by about 0.1–0.3%. Consequently, we give a conservative estimate of the bias on the polarization fraction percentiles (and therefore on p_max) by considering the difference $\overline{h(p_i)}-h(p)$. A rough debiasing of the data percentiles by this quantity is achieved by subtracting this value from h(p). For instance, the estimated bias for the 99.9th percentile at 80′ resolution in the fiducial offset case is about 0.66%. Subtracting this from the data percentile, we obtain a debiased value of 22.00%.

Finally, we emphasize that the truly dominant source of uncertainty in the determination of characteristic values of the p distribution is the offset in I. It is larger than the statistical uncertainty, which is of order 0.01–0.10%, or the impact of the residual systematics that has been estimated in Appendix C to be typically 0.5%.

Performing the same debiasing for the low and high offset values, and gathering these results for the 99.9th percentile, we obtain a debiased value of $22.0_{-1.4}^{+3.5}\pm 0.1\%$ for the maximum dust polarization fraction observed at 80′ resolution and 353 GHz over the full sky, where the first uncertainty relates to the systematic effect of the total intensity offset and the effects of residual systematics, and the second covers the statistical uncertainty estimated from the 1000 Monte Carlo realizations.

For completeness, Table 1 also gives the same percentiles for the maps smoothed to 160′ resolution, showing a further reduction of the bias $\overline{h(p_i)}-h(p)$. In that case, we find that the maximum dust polarization fraction observed is $21.4{}_{-1.2}^{+2.2}\pm 0.1\%$. This value of p_max and the debiased value at 80′ agree quite well. This shows that smoothing has little effect on the polarization fraction. Of course, the amount of smoothing applied should not be excessive, because of the potential impact of beam depolarization at higher FWHM. In Appendix F.8, we quantify the effect of smoothing on p and p_max in the framework of the analytical model presented in Sect. 4.3. It is found that smoothing from one resolution to another leads to a decrease in p² by an amount that is statistically independent of the value of the polarization fraction. Considering p itself, this means that the effect of smoothing is very small if p is large, e.g., p ≈ p_max (Appendix F.9). We conclude that our derivation of p_max is not so much affected by the resolution and much more so by the offset in I.

These results are consistent with the finding of Planck Collaboration Int. XIX (2015) that p_max >  19.8% at 60′ resolution over a smaller fraction of the sky. We have also checked that they are not significantly affected when selecting only those pixels on the sky for which the S/N in polarization is p/σ_p >  3.

As was pointed out in Planck Collaboration Int. XIX (2015) and Planck Collaboration Int. XX (2015), the level of observed polarization fractions is strongly dependent on the angle Γ of the mean magnetic field with respect to the plane of the sky (see Appendix F and Fig. F.2). The distribution function of p must depend on this mean orientation of the Galactic magnetic field. Compared to what would be obtained for a mean field that is everywhere in the plane of the sky, the distribution should be more peaked towards lower values, as we do observe, but the value of p_max might still be high, reflecting those parts of the sky with a favourable orientation, i.e., in the plane of the sky. Although the estimate of p_max based on percentiles would be impacted, such an analysis (requiring a model of the large-scale GMF) is beyond the scope of this paper.

4.1.2. Polarization angle

Figure 6 shows the distribution function for the polarization angle ψ, for which the value of the total intensity offset is unimportant. The comparison between the histogram for the GNILC map of ψ and the mean histogram over the Monte Carlo realizations shows that there is virtually no noise bias. The histograms peak around 0°, which corresponds to an orientation of the GMF parallel to the Galactic plane. Quantitatively, over the 1000 Monte Carlo samples, the ensemble average of the mean polarization angle is −0${}_{.}^{\circ }$64 ± 0${}_{.}^{\circ }$03. This value is compatible with the earlier measurement in Planck Collaboration Int. XIX (2015, see their Fig. 3).

4.1.3. Polarization angle dispersion function

Finally, the distribution function of 𝒮 is shown in Fig. 7. Results for the case of a 160′ FWHM and lag δ = 80′ are shown in green, and for the case of a 80′ FWHM and lag δ = 40′ in blue. As above for p and for ψ, the solid lines are for the GNILC maps, the dashed lines are the Monte Carlo means, and the coloured regions show the span of histograms for the 1000 Monte Carlo realizations.

It is interesting that these distributions have a tail passing through 52°, the value of 𝒮 for randomly oriented polarization. As noted by Alina et al. (2016), if an orientation distribution produces a value of 𝒮 that is somewhat lower (higher) than this, then the addition of noise tends to make 𝒮 larger (smaller), towards 52°. This tail in the full DF in Fig. 7 is strongly associated with regions where p is small and more susceptible to the influence of noise, as is apparent in Fig. 8, which shows the distribution function of 𝒮 for different ranges in polarization fraction (p <  1%, 1%< p <  5%, and p >  5%) for the GNILC data at 160′ resolution and with a lag δ = 80′. The large values of 𝒮 are also associated with large values of the scatter ${\sigma }_{\mathcal{S}}^{\mathrm{MC}}$, as shown by the widening of the envelope at high values of 𝒮 in Fig. 7. The width of the envelope at 160′ resolution is compatible with the largest values found in the map of ${\sigma }_{\mathcal{S}}^{\mathrm{MC}}$ (Fig. 4).

Figure 7 shows that, for the case of an 80′ FWHM and lag δ = 40′, at large values of 𝒮 the mean DF of the Monte Carlo realizations is clearly biased with respect to the distribution function of the data, which does not even fit within the region spanned by the 1000 Monte Carlo simulations. On the other hand, for 160′ FWHM and lag δ = 80′, the bias is much less apparent and so we focus on the results for this case. Despite the long tail at large 𝒮, most of the points in this tail have low occurrence rates, underlining the regularity of the polarization angle on large scales. At 160′ resolution and a lag of δ = 80′, the distribution of values in the 𝒮 map for the data peaks at 1${}_{.}^{\circ }$7, with mean and median values of 7${}_{.}^{\circ }$6 and 4${}_{.}^{\circ }$6, respectively. The same characteristic values over the 1000 Monte Carlo simulations are, respectively, 1${}_{.}^{\circ }$9 ± 0${}_{.}^{\circ }$6, 8${}_{.}^{\circ }$29 ± 0${}_{.}^{\circ }$01, and 5${}_{.}^{\circ }$12 ± 0${}_{.}^{\circ }$01. Using the 99th percentile, most of the points in the data have 𝒮 ≤ 43${}_{.}^{\circ }$6, while the Monte Carlo simulations give an estimate of 45${}_{.}^{\circ }$3 ± 0${}_{.}^{\circ }$2. We give these values for reference in the future, for instance in work comparing Planck data with MHD simulations and analytical models.

We stress that while the smoothing to 160′ is warranted here for studies including the high-latitude sky, this requirement for smoothing should not be generalized. Indeed, when the analysis is restricted to the approximately 42% of the sky considered in Planck Collaboration Int. XIX (2015), we find that no such bias exists when working at 80′ FWHM and lag δ = 40′. Incidentally, this confirms the results shown in Planck Collaboration Int. XIX (2015) at 60′ resolution and δ = 30′.

4.2. Two-dimensional distribution functions

In this section we investigate the 2-dimensional joint distribution functions of polarization fraction p and another variable. Therefore, instead of simply presenting a scatter plot, we display a 2-dimensional histogram made by binning in the two dimensions and encoding the number in each bin by colour.

4.2.1. Polarization fraction versus total gas column density

In Fig. 9 we display the 2-dimensional histogram of p and total gas column density N_H, using the GNILC data, at 353 GHz and uniform 80′ resolution, with the fiducial total intensity offset, over the full sky. We discuss the determination of N_H at the beginning of Sect. 5. No cut in either S/N or Galactic latitude has been performed here. The colour scale encodes the logarithm of counts in each bin, while the black curves show the 5th, 95th, and 99th percentiles of the p distribution in each N_H bin, as well as the median polarization fraction in each N_H bin.

To explore the sensitivity of this distribution and characteristic curves to statistical noise, we use the Monte Carlo approach described above. We first compute the 2-dimensional distribution function of p and N_H for each of the 1000 simulations, along with the curves giving the median and the 5th, 95th, and 99th percentiles of p within each bin of N_H. We then compute the average curve for each of these four quantities, as well as their dispersions within each N_H bin.

We find that these exhibit small statistical dispersions, but that towards the most diffuse lines of sight (N_H <  10²⁰ cm⁻²), the maximum polarization fractions (measured for instance by the 99th percentile curve) for the Monte Carlo simulations are slightly higher than the corresponding values from the data. As expected, this bias is in the same sense as discussed in Sect. 4.1.1 for the distribution function of p. Recall that this is for 80′ resolution; when working at 160′ resolution this bias disappears.

The joint (N_H,p) distribution has qualitatively the same behaviour as that found over a smaller fraction of the sky in Planck Collaboration Int. XIX (2015): a large scatter of p towards diffuse lines of sight and a decrease in the maximum p as N_H increases.

For completeness, we show in Appendix E the effect of the total intensity offset. It is negligible at the high intensity end, where the histograms are similar whether we use the fiducial, high, or low offset values. At the low intensity end, on the other hand, the effect of the offset is more marked. There is a significant increase in characteristic values (highest percentiles) of p for decreasing N_H when taking the low offset, and conversely a marked decrease in the maximum p when taking the high offset.

One might wonder if it would be possible to constrain the offset by assuming that p should reflect dust properties at low column densities, and therefore that the offset should be such that the maximum p is approximately constant at low N_H. In this respect, the fiducial offset seems more adequate than either the high or low cases, as can be seen by comparing Fig. 9 with Fig. E.1.

The sharp downturn of the maximum polarization fraction observed near N_H ≈ 10²² cm⁻² corresponds to the strong depolarization occurring on lines of sight that probe high column density structures that are not resolved at 80′.

4.2.2. Polarization angle dispersion versus polarization fraction

In Planck Collaboration Int. XIX (2015), we discovered an inverse relationship between the polarization fraction p and the polarization angle dispersion function 𝒮, working with data over approximately 42% of the sky, at a resolution of 60′ and a lag of δ = 30′. We have verified quantitatively on the same sky region and using the same methodology that the same inverse relationship holds with the Planck 2018 data release; the maps of polarization are very similar where the S/N is high, as expected. In this limited sky region, we also find that the results are only slightly dependent on the adopted Galactic offset.

Extending to the full sky at 160′ resolution and a lag of δ = 80′, we present the 2-dimensional histogram of the joint distribution function of 𝒮 and p in Fig. 10. The data clearly show that the inverse relationship seen at low and intermediate Galactic latitudes in Planck Collaboration Int. XIX (2015) persists in the high-latitude sky. In Fig. 10 we also display the running mean of 𝒮 in each bin of p for the data¹³. We show in the next section that simple analytical models suggest that the relationship is indeed ⟨𝒮⟩_p ∝ 1/p. Such a trend is shown in Fig. 10 as the dashed white line.

4.3. Relationship to models

All of the properties discussed so far, namely the distribution functions of p, ψ, and 𝒮, the decrease in the maximum polarization fraction with increasing column density, and the inverse relationship between 𝒮 and p, are consistent with the analysis first presented in Planck Collaboration Int. XIX (2015). Subsequently, phenomenological models of the polarized sky incorporating Gaussian fluctuations of the magnetic field have been developed (Planck Collaboration Int. XLIV 2016; Ghosh et al. 2017; Vansyngel et al. 2017; Levrier et al. 2018). Interestingly, although these models were built to reproduce some 1- and 2-point statistics of polarized emission maps, they were not tailored to reproduce the inverse relationship between 𝒮 and p that is evident in the Planck data, and yet they are able to do so very readily and robustly. A similar inverse relationship between 𝒮 and p was also observed in synthetic polarization maps built from numerical simulations of MHD turbulence (Planck Collaboration Int. XX 2015).

In Appendix F we present an analytical derivation of this property, based on the most basic version of these phenomenological models. In that framework, the emission is assumed to arise from a small number N of layers, each emitting a fraction 1/N of the total intensity, and harbouring a magnetic field that is the sum of a uniform component and a turbulent Gaussian component. The main parameters of the model, besides N, are the intrinsic polarization fraction p₀¹⁴, the level of the turbulent magnetic field f_M relative to the magnitude of the uniform component, and the spectral index α_M of the spatial distribution of this turbulent component.

In Fig. 11 we show the 2-dimensional distribution function of 𝒮 and p at 160′ resolution, using a lag δ = 80′, for a polarized sky built from such a Gaussian phenomenological model, with α_M = −2.5, f_M = 0.9, N = 7, and p₀ = 26%. This choice of parameters, within the range of good fits in Planck Collaboration Int. XLIV (2016), leads to the mean analytical relation (see Appendix F.9) ⟨𝒮⟩_p = 0${}_{.}^{\circ }$34/p, which is very close to a fit to the observational trend, ⟨𝒮⟩_p = 0${}_{.}^{\circ }$31/p, overplotted in Figs. 10 and 11. We show in Appendix F.9 that this relationship depends weakly on the resolution ω according to

$$\begin{array}{c}\hfill {\langle \mathcal{S}\rangle }_{\mathrm{p}}=\frac{0.{}^{°}31}{p}{\left(\frac{\omega }{{160}^{\prime }}\right)}^{0.18}.\end{array}$$(10)

Because changes of dust properties or dust alignment are not included in these phenomenological models nor in the synthetic observations from MHD simulations, we conclude that the inverse relationship between 𝒮 and p is a generic statistical property that results primarily from the topology of the magnetic field.

We also note that neither the phenomenological model of Planck Collaboration Int. XLIV (2016), nor the MHD simulation in Planck Collaboration Int. XX (2015), account for the 3D structure of the ordered (mean) component of the GMF on large scales. The imprint of this ordered component on the dust polarization can be readily identified in the map of the dust polarization angle (Fig. 2). It also impacts the polarization fraction map on large angular scales and thereby the dependency of p on Galactic latitude. For synchrotron polarization, this has been quantified by Page et al. (2007) and Miville-Deschênes et al. (2008) using Galaxy-wide models of the GMF. As discussed in Alves et al. (2018), a comprehensive model of dust polarization would also need to take into account the structure of the GMF on the scale of the Local Bubble (100–200 pc).

5. Insight from interrelationships and Galactic context

Further insight into statistical measures of the polarization can be gained not only by considering them in relation to one another, but also by studying how they jointly vary with other physical parameters such as dust temperature T_d or column density, and how these relationships evolve from the diffuse ISM to molecular clouds.

An important parameter in this study is the total amount of dust along the line of sight, or dust column density, which is usually quantified by the dust optical depth τ (at 353 GHz). Because dust emission is optically thin at this frequency, this relates the modified blackbody (MBB) model of the emission to the total intensity via

$$\begin{array}{c}\hfill I_{\nu }=\tau B_{\nu }(T_{\mathrm{d}}){\left(\frac{\nu }{353\mathrm{GHz}}\right)}^{\beta },\end{array}$$(11)

where β is the observed dust spectral index (Planck Collaboration XI 2014). It is also common to rescale from dust optical depth to entirely different units like colour excess in the optical, E(B − V)_τ¹⁵, or total column density of hydrogen N_H. The calibrations of such rescalings are uncertain and possibly dependent on the environment. This is not important for our results below and we use the MBB parameters τ and T_d from Planck Collaboration Int. XLVIII (2016), the calibration from Planck Collaboration XI (2014) at 353 GHz,

$$\begin{array}{c}\hfill E{(B-V)}_{\tau }=(1.49\pm 0.03)\times {10}^4\tau ,\end{array}$$(12)

and the relation N_H = 5.87 × 10²¹ cm⁻² × E(B − V)_τ (Bohlin et al. 1978; Rachford et al. 2009) to estimate N_H. It is preferable to use τ converted to N_H rather than an estimate of the gas column density derived from H I because of the presence of dust in the WIM that is sampled by all of our polarized and unpolarized observables.

We note that in this section we use not only the GNILC maps at 80′ and 160′ resolution, but also the alternative Stokes maps (Sect. 2.4) at finer resolutions of 40′, 20′, and 10′.

5.1. Origin of the observed variations of the polarization fraction p

The mutual correlations between p, 𝒮, and the column density N_H were studied in detail for the particular case of the Vela C molecular cloud by Fissel et al. (2016) using BLASTpol data. From the present Planck 2018 data, Fig. 12 shows how these correlations appear for the more diffuse ISM (4 × 10¹⁹ cm⁻² <  N_H <  10²² cm⁻²) over the whole sky, excluding only the latitudes close to the Galactic plane (|b| < 5°). Significant variations about the trend of p with N_H prevent modelling it by a simple relationship. For N_H <  5 × 10²⁰ cm⁻², the mean value is compatible with a constant, then decreases over the range 0.5–2 × 10²¹ cm⁻², and eventually becomes rather flat again. Colouring¹⁶ the data with 𝒮 (on a logarithmic scale) we see from the stratification of the data in Fig. 12 that there is a gradient of 𝒮 mainly perpendicular to the observed trend of p with respect to N_H. This analysis indicates that the decreases in p with 𝒮 and with N_H are mostly independent of each other.

Figure 13 shows how the variations of p and 𝒮 at a given column density are related to the dust temperature T_d. Dust tends to be systematically cooler when p is high and warmer when p is low (top panel). This is observed at all but the highest column densities (N_H >  8 × 10²¹ cm⁻²). On the other hand, the opposite is seen in 𝒮 (bottom panel). The mirror symmetry between the two panels of Fig. 13 shows convincingly that there is in fact no physical relation between the polarization fraction p and the dust temperature T_d in the diffuse ISM. Even if it seems that, at any column density, high p corresponds to colder dust and low p to warmer dust, the bottom panel demonstrates that the value of p is actually driven by 𝒮, i.e., by the magnetic field structure and the depolarization produced by its variations along the line of sight and within the beam.

5.2. Exploring beyond first-order trends using 𝒮  ×  p

In Sect. 4.3, we concluded that the inverse relationship between 𝒮 and p is a generic statistical property that results from the topology of the magnetic field alone, and that a trend 𝒮 ∝ 1/p, close to that observed, is expected on the basis of simple analytical models (Appendix F). It is therefore interesting to explore beyond this underlying cause for the inverse relationship, in search of evidence for the impact of other physical factors, such as dust alignment efficiency, elongation, and composition. For this we can use the product 𝒮 × p, which removes the impact of the magnetic field structure statistically. This does not mean that the product depends only on properties of the dust, e.g., the maximum polarization fraction p_max. As explained in Appendix F, the product 𝒮 × p also depends on the length over which dust structures are probed along the line of sight, as well as on the ratio of the turbulent to ordered magnetic field. Nevertheless, it is interesting to try this approach, as also emphasized by the mirror symmetry seen between the two panels of Fig. 13.

Accordingly, Fig. 14 compares the variations of not only p and 𝒮, but also 𝒮 × p with N_H, Galactic latitude b, and Galactic longitude l. It should be noted that throughout this entire analysis lines of sight close to the Galactic plane (|b| < 5°) are excluded.

As expected, the product 𝒮 × p has smaller variations with N_H, b, and l than exhibited by p and 𝒮 separately, and the decrease in 𝒮 × p with N_H is systematic, without significant departures.

Away from the Galactic plane, the dependence of 𝒮 × p on b is less pronounced than it is for p and also more symmetric between positive and negative latitudes. The strong dependence of p on b that can be attributed to the systematic change in the orientation of the mean magnetic field with respect to the line of sight is mitigated in 𝒮 × p, confirming our interpretation. However, there are still small variations of 𝒮 × p over a large spatial scale that remain to be interpreted. Towards the Galactic plane there is a pronounced dip, that is probably due to the accumulation of variously polarized structures along the line of sight at these low latitudes (Jones et al. 1992). This dependence on the latitude will be further discussed in Sect. 5.4.

As with the dependence on N_H, the variations of the product 𝒮 × p with Galactic longitude l are much less pronounced (of the order of 30%) than those of p and 𝒮 independently (a factor 3 or so in each case).

5.3. Dedicated study for six molecular regions in the Gould Belt

The radiative torques theory (RAT; Lazarian & Hoang 2007; Hoang & Lazarian 2016) makes strong predictions for the dependence of the dust alignment on local physical conditions, namely the intensity and anisotropy of the radiation field, and the angle between the magnetic field and the anisotropic radiation field. Dense regions, screened from the interstellar radiation field and possibly with embedded sources, should be promising regions in which to identify evidence for RATs (Vaillancourt & Andersson 2015; Wang et al. 2017).

To probe this possibility, we have selected a set of six 12° × 12° molecular regions in the Gould Belt (Dame et al. 2001). These are listed in Table 2. All but one (Aquila Rift) were already studied using Planck 2013 data in Planck Collaboration Int. XX (2015). The higher S/N in these bright, high column density regions enables an analysis at a higher resolution (40′, N_side = 256), and the uncertainty on the offset in total intensity I can be safely ignored. For this study, we therefore make use of the ASMs (Sect. 2.4) with the fiducial total intensity offset.

Figure 15 (top panel) shows that the variation of p with N_H is very diverse in these molecular clouds, as was already observed in Planck Collaboration Int. XX (2015). For some (e.g., Aquila Rift and Chamaeleon) p is fairly high in the more diffuse envelope but progressively decreases towards denser parts. Others (e.g., Polaris and Orion) have low p at all column densities, while for one (Ophiuchus) p increases, then decreases. In each region, the corresponding variations of 𝒮 with N_H (Fig. 15, middle panel) are clearly inversely related to those of p, so that the product is by contrast almost constant and uniform across the sample of clouds, as can be seen in Fig. 15 (bottom panel).

Figure 16 (top panel) directly shows this inverse trend between 𝒮 and p in the selected molecular clouds. The various curves show the mean 𝒮 in each bin of p, all bins containing the same number of pixels. Also plotted is the average curve over all the different regions. Despite differences in mean column densities, regions as different as Polaris and Orion all fall on the same correlation line. This figure demonstrates that most of the variations of p with N_H in Fig. 15 can be attributed to variations of 𝒮 alone, and eventually to the variation of the magnetic field structure along the line of sight and in the plane of the sky.

We note that the inverse relationship between 𝒮 and p is the same as the one found in Sect. 4.2.2 over the full sky (Fig. 10), which is dominated by the high-latitude diffuse ISM, once the difference in resolution, and therefore in the lag δ used to compute 𝒮, is accounted for in the framework of our analytical model (see Appendix F.7). Indeed, for this analysis towards selected molecular regions, we work at a finer 40′ FWHM resolution, instead of 160′, and numerical results of the model show that the product 𝒮 × p scales as FWHM^0.18 (Eq. (10)). The prediction of this model is shown as the dashed lines in Figs. 15 and 16.

The bottom panel of Fig. 16 shows the result of the same procedure applied to the phenomenological model described in Sect. 4.3 and smoothed to the same resolution of 40′. Our model, which was able to reproduce the trend 𝒮 ∝ 1/p observed at large scale (Figs. 10 and 11), can also reproduce it at the smaller scales probed here, in regions of 12° × 12°. The downward shift of the correlation observed in the data (compare Fig. 16 with Figs. 10 and 11), that is due to the change in resolution and is already integrated in our expression for 𝒮 × p, is also observed in the simulation.

As quantified in Appendix F within the framework of the phenomenological model of Planck Collaboration Int. XLIV (2016), the mean value of the 𝒮 × p product depends on f_M, the ratio of the dispersion of the turbulent component of the field to the mean field strength (see Eq. (F.53), where f_m scales linearly with f_M). Thus, the alignment of the data lines in the top panel of Fig. 16 is a remarkable result, which suggests that the strength of the turbulent component of the magnetic field, relative to the mean field strength, is comparable among Gould Belt clouds, and between clouds and the diffuse ISM, despite differences in the local star-formation rate. This interpretation is illustrated by the correspondence between the top and bottom (data and model) plots of Fig. 16.

In the cold neutral medium, the magnetic and turbulent kinetic energies are known from H I Zeeman observations to be in approximate equipartition (Heiles & Troland 2005). Our analysis of the Planck data suggests that this energy equipartition also applies to the Gould Belt clouds. This result is consistent with the much earlier results derived from the modelling of stellar polarization data by Myers & Goodman (1991) and Jones et al. (1992).

To test the RAT theory, we need to estimate the relative intensity of the radiation field, G₀, in these regions and then look for a possible correlation between this value and the polarization fraction, once the latter is adjusted for the variations related to 𝒮, i.e., look for a correlation between G₀ and 𝒮 × p. To this end, we use an estimator G_ℛ of the radiation field intensity (Guillet et al. 2018; Fanciullo et al. 2015) that is based on the assumption of thermal equilibrium for large dust grains, which dominate the emission at this frequency. Under this hypothesis, the dust radiance ℛ, which is the integrated intensity of the dust emission (Planck Collaboration XI 2014), is balanced by the heating of dust by absorption of the ambient radiation field. The relative intensity G₀ of this ambient field is therefore estimated using the radiance per unit visual extinction A_V, and in practice, the estimate G_ℛ is computed through

$$\begin{array}{c}\hfill G_{\mathcal{R}}=\frac{E{(B-V)}_{\mathcal{R}}}{E{(B-V)}_{\mathrm{Green}}},\end{array}$$(13)

where E(B − V)_ℛ stands for the dust radiance converted to a colour excess using a correlation with extinction to quasars (Planck Collaboration XI 2014), while E(B − V)_Green is from a colour excess map (Green et al. 2018) based on Pan-STARRS1 (PS1).

Figure 17 shows the correlation between this estimate G_ℛ of the radiation field intensity and the dust temperature T_d from the MBB fit, in the molecular regions selected, with data points coloured according to the dust optical depth τ converted to a column density N_H. This is the equivalent of Fig. 2 in Guillet et al. (2018), limited to the Gould Belt regions selected. We see that the correlation is quite tight and follows a scaling $T_{\mathrm{d}}\approx 18.5{G_{\mathcal{R}}}^{1/5.6}$, corresponding to an average temperature of 18.5 K for a standard radiation field G₀ = 1 and a spectral index β = 1.6. We note that red points at low dust temperatures do not follow this trend perfectly because the reddening map of Green et al. (2018) tends to underestimate the true column density at high optical depths, so that G_ℛ is overestimated. We also note that at high optical depths the spectral shape of the ambient radiation field is altered by the frequency dependence of the extinction, which would also impact the amount of energy absorbed and thus the dust temperature. Overall, over a wide range of G_ℛ this plot demonstrates that a change in the dust temperature is a good indicator of a change in the ambient radiation field intensity G₀ in these molecular regions.

According to the RAT theory, we would expect grains in a more intense radiation field to be better aligned, and therefore that 𝒮 × p would tend to increase with increasing T_d. However, Fig. 18 does not show any correlation between the polarization fraction and the dust temperature, and the product 𝒮 × p, which we use as a proxy for dust alignment, does not show any trend with T_d either. This analysis for molecular clouds at a resolution of 40′ confirms our conclusion drawn in Sect. 5.1 for the diffuse ISM that there is no strong indication of a link between the polarization fraction p and the dust temperature T_d.

5.4. Multi-resolution view of the variations of p, 𝒮, and 𝒮  ×  p across the ISM

The discussion of the variations of p, 𝒮, and 𝒮 × p in the diffuse ISM at 160′ (Sect. 5.2) and in molecular clouds at 40′ (Sect. 5.3) suggests a multi-resolution exploration of these trends. In Fig. 19, we present such a view by compiling mean trends at all resolutions, from 10′ to 160′. The impact of noise bias at low column densities in p, and even more so in 𝒮, is clearly seen as a rising deviation (dotted line) from the bundle of curves that otherwise roughly match – except for a global shift – at higher column densities, despite the fact that both p and 𝒮 have been debiased. All debiasing methods indeed fail when the S/N becomes lower than 1. This occurs below a threshold in N_H that is different for p and 𝒮 and increases with decreasing FWHM. Our debiasing methods are known to be statistically efficient when the S/N is higher than about 3 (Simmons & Stewart 1985; Montier et al. 2015a).

Figure 20 presents the evolution of the mean S/N in p and 𝒮 as a function of the column density for various resolutions of the map (still considering only |b| > 5°). Obviously, the S/N tends to increase with N_H at a given resolution, and to decrease with increasing resolution at a given N_H. With these figures we can estimate, at each resolution, the column density threshold above which the mean S/N for p and 𝒮 is larger than 3, i.e., above which debiased values of p and 𝒮 are robust. The S/N is smaller for 𝒮 than it is for p, and we note the large scatter in S/N for any given bin in N_H. Users of Planck data in polarization should take into account these thresholds in column density to estimate the reliability of the debiasing.

In Fig. 19, the data points below these thresholds are excluded from our analysis (dotted curves). The spread of values for p and 𝒮 × p at low column densities, induced by the uncertainty on the offset in total intensity I, is indicated by dashed lines for 160′. For reference, the density of points at a resolution of 10′ is plotted as a background. It is only plotted for N_H >  4 × 10²¹ cm⁻², which is the column density threshold for 𝒮 at this resolution.

Inspection of Fig. 19 shows that there always exists a range in column density where the curves at two consecutive resolutions are parallel to each other, i.e., they probe the same variations. It is therefore possible to obtain a unique, smooth and continuous trend for each quantity as a function of N_H through a renormalization of the profiles at each resolution, leading to Fig. 21. We proceed as follows, for p, 𝒮, and 𝒮 × p. At each resolution ω (expressed in arcmin), using Fig. 20 we determine the minimal column density ${N_{\mathrm{H}}}^{\min }(\omega )$ above which 𝒮 is correctly debiased, and define a reference interval $\mathcal{I}(\omega )=[{N_{\mathrm{H}}}^{\min }(\omega ),10\times {N_{\mathrm{H}}}^{\min }(\omega )]$ (indicated by the horizontal colour bars in Fig. 21). For two consecutive resolutions ω and ω/2, we compute the mean values of p, 𝒮, and 𝒮 × p at both resolutions on the common interval ℐ(ω)∩ℐ(ω/2) and then the ratio of these two values, r_ω, ω/2. Finally, we compute the factor by which each curve at resolution ω from Fig. 19 must be multiplied to be normalized to the curve at the coarsest resolution, ω_max = 160′, i.e., r(ω) = r_{ωmax, ωmax/2} × ... × r_4ω, 2ω × r_2ω, ω. This renormalization removes the depolarization induced in p by the smoothing of the data, as well as the change of the lag δ with the resolution, as far as 𝒮 in concerned.¹⁷

The mirrored similarity of each detailed variation in the logarithmic representation of p and 𝒮 in Fig. 21 clearly emphasizes the inverse relationship between these two quantities. In our multi-resolution normalized representation of the variations of p with the column density, the mean value of p decreases by a factor 3−4 from the lowest column densities at high latitudes and a resolution of 160′, to the highest column densities probed here at a resolution of 10′. This strong decrease is almost entirely mirrored as an increase in 𝒮 by the same factor, demonstrating again that the decrease in p with the column density is mainly a consequence of the depolarization by the structure of the magnetic field.

The residual trend in 𝒮 × p is small, about a 25% decrease with column density from 10²⁰ cm⁻² to 2 × 10²² cm⁻². For the case |b| > 10°, the profile of 𝒮 × p over this same range of N_H is quite flat. For the case including |b| > 2°, 𝒮 × p decreases systematically with N_H. In our phenomenological model, a systematic decrease in 𝒮 × p is expected at low Galactic latitudes from an increased number of independent layers N along the line of sight (see Eq. (F.45)), related to the increased dust opacity and/or length probed (Jones et al. 1992), that is not compensated by the inverse effect of increased 𝒮 due to an increased distance to the observed dust material (recall that 𝒮 depends on the physical scale probed in the cloud, therefore on its distance). There remains therefore little room for a systematic variation of grain alignment for column densities up to 2 × 10²² cm⁻².

At slightly higher column densities (N_H >  3 × 10²² cm⁻²), we observe a decrease in p, together with an increase in 𝒮 and 𝒮 × p. Such a combination cannot be explained by a decrease in grain alignment, which would not affect 𝒮. These lines of sight, part of the Orion and Ophiuchus regions, are situated at intermediate latitudes (10° <  |b| < 20°) and probably do not suffer from depolarization by the Galactic background, unlike other lines of sight at lower column densities situated at lower latitude. As can be seen from the density of points in Fig. 19, this departure from the mean trend has a low statistical significance, which prevents us from commenting further.

To conclude, most variations of the polarization fraction p with N_H are inversely related to those of 𝒮, a tracer of the depolarization by the turbulent magnetic field. The multi-resolution study of the variation of 𝒮 × p with N_H does not indicate any systematic variation of the grain alignment efficiency beyond around 25%, up to a column density of 2 × 10²² cm⁻².

5.5. Grain alignment efficiency in the ISM

In this section, we discuss the impact of our results on the question of grain alignment.

Since the pioneering work of Myers & Goodman (1991) and Jones et al. (1992), it has been clear that the structure of the magnetic field along the line of sight plays a major role in the build-up of polarization observables. Nevertheless, the decrease in the polarization fraction with increasing column density is widely considered as evidence for a systematic decrease in the degree of grain alignment efficiency with increasing exinction (Whittet et al. 2008; Cashman & Clemens 2014; Alves et al. 2014).

In this paper, we have shown that most (if not all) variations observed in the polarization fraction p are mirrored in the dispersion of polarization angles, 𝒮, a quantity that is independent of the grain alignment efficiency and of dust optical properties. Quantitatively, the near constancy of 𝒮 × p with increasing column density indicates that the variations of the polarization fraction are dominated by the structure of the magnetic field, not only in the diffuse ISM, but also in molecular clouds, at least up to a column density of N_H ≈ 2 × 10²² cm⁻². The decrease in grain alignment efficiency with column density cannot exceed about 25%, from the most diffuse ISM up to this same column density, N_H ≈ 2 × 10²² cm⁻². These results are significant constraints for theories of grain alignment.

Dust alignment in the ISM is widely thought to be associated with radiative torques (RATs). As mentioned, in the classical framework of RATs (Lazarian & Hoang 2007), the grain alignment efficiency depends on the radiation field intensity and on the angle between the radiation field anisotropy and the magnetic field. During the last decade, there have been a few claims for evidence of such effects (e.g., Andersson & Potter 2010; Vaillancourt & Andersson 2015; Andersson et al. 2015). Analysing Planck full-sky data, we could not find, either in the diffuse ISM or in molecular clouds, any signature in polarization observables that could point to a significant variation of grain alignment related to a variation in the grain temperature.

However, the low resolution of Planck data (5′), combined with the smoothing of the maps necessary to guarantee that p and 𝒮 are not biased by noise (160′ in the high latitude diffuse ISM, 40′ in molecular clouds), does not allow us to probe the same physical conditions as, for example, NIR polarimetry through dense clouds (Jones et al. 2015). A detailed analysis of Planck polarization without the additional smoothing, and hence at higher resolution (7′) and higher column densities, will be pursued in a future paper dedicated to cold cores.

6. Comparison with starlight polarization at high Galactic latitudes

In this section, we correlate Planck polarization with starlight polarization in the diffuse ISM at high Galactic latitudes, and derive new constraints on dust models. Following the approach in Planck Collaboration Int. XXI (2015) for translucent lines of sight (0.15 <  E(B − V)_τ <  0.8) at low Galactic latitudes (|b| < 30°), we can now derive emission-to-extinction polarization ratios for the diffuse high-latitude Galactic ISM (E(B − V)_τ <  0.15, corresponding to column densities N_H <  10²¹ cm⁻²). The ratios are

$$\begin{array}{cc}& R_{P/p}=\frac{P}{p_V};\hfill \end{array}$$(14)

$$\begin{array}{cc}& R_{\mathrm{S}/\mathrm{V}}=\frac{P/I}{p_V/{\tau }_V}·\hfill \end{array}$$(15)

Here P and I are what Planck has measured in the submillimetre. In the optical V band, p_V is the degree of polarization for a star to which the optical depth is τ_V. The latter is estimated from the colour excess of the star, E(B − V)^⋆, using A_V = R_V × E(B − V)^⋆ with the ratio of total to selective extinction R_V = 3.1 (Fitzpatrick & Massa 2007), and τ_V = A_V/1.086.

These polarization ratios quantify the amount of polarized emission per unit polarized extinction. Because they measure the effects of the same grains at different wavelengths, they remove the first-order dependencies of the polarization observables on the magnetic field structure and grain alignment efficiency (Planck Collaboration Int. XXI 2015). As such, they directly provide observational constraints on the optical properties of the dust population that is aligned with the magnetic field, and strongly constrain dust models (Guillet et al. 2018).

The second of these ratios, R_S/V, being inversely proportional to I, is sensitive to the total intensity offset. We comment on the derived values of R_S/V for the low, fiducial, and high values of the offset in Sects. 6.4 and 6.5.

We are interested in examining how the amount of submillimetre polarization is related to the amount of optical polarization from the same dust. If the dust probed in the submillimetre and the optical is the same, then the polarization orientations should be orthogonal. As discussed in Appendix G.4, this is the case for the lines of sight used in our analysis of the ratios. This is quantified by the polarization angle difference Δψ_S/V which takes into account the 90° difference (see Eq. (G.6)).

In Appendix G we describe in detail our estimates of the noise and systematic uncertainties that affect the submillimetre and optical data as used in this new analysis. We highlight the relevant results in the discussion below.

6.1. Estimates for starlight reddening

To enable appropriate comparison of polarization properties in the optical and in the submillimetre on lines of sight to stars, it is necessary to obtain an estimate of the reddening to the star. To this end, we obtain the distance D^⋆ to each star by extracting the parallax θ^⋆ from the Gaia DR2 release (Gaia Collaboration 2016, 2018) or from the polarization catalogues when Gaia data are not available. Then we derive an estimate of the reddening to the star, E(B − V)^⋆, by interpolating the PS1-based 3-dimensional reddening data cube (Green et al. 2018) at distance D^⋆. This data cube is composed of 31 maps, each representing a range in distance modulus. To limit the impact of noise in our analysis, the 31 maps were separately smoothed to a resolution of 40′ and downgraded to N_side = 256. We also converted the PS1-based reddening to the Johnson E(B − V) scale using the relation E(B − V)_Johnson = E(B − V)_PS1/1.0735 (Table 6 of Schlafly & Finkbeiner 2011). The Johnson scale is used hereafter without explicit subscripting. From the same maps we also obtain the total reddening along the line of sight, E(B − V)^∞. Uncertainties related to the reddening maps are estimated in Appendix G.3.

6.2. Polarization data

For this analysis, we use the alternative Planck 353-GHz I, Q, and U maps from the Planck 2018 data release (ASMs, Sect. 2.4), smoothed to a resolution of 40′ to limit the noise in Q and U.

From a series of optical polarization catalogues of high-latitude stars (Berdyugin et al. 2001, 2004, 2014; Berdyugin & Teerikorpi 2001, 2002), we extract data for 2461 lines of sight to stars with measured degree of polarization p_V, uncertainty σ_pV, and polarization angle ψ_V (in the IAU convention, consistently with our definition of ψ for the Planck data). These catalogues cover the northern Galactic hemisphere at high latitudes (2092 stars with b >  30°), and part of the southern cap (369 stars with b <  −59°). After removing 70 stars falling outside the region covered by PS1 (mainly in the southern hemisphere) and 3 stars without a distance estimate, there remain 2388 stars for which we have both reddening estimates and optical polarization data.

As with the Planck submillimetre data, we use the MAS estimator (Plaszczynski et al. 2014) to debias the degree of polarization p_V in the optical. Using these values of p_V and ψ_V, we then define Stokes parameters in extinction, q_V and u_V, in the same HEALPix convention as our Planck data:

$$\begin{array}{c}\hfill q_V=p_{V}cos(2{\psi }_V);\end{array}$$(16)

$$\begin{array}{c}\hfill u_V=-p_{V}sin(2{\psi }_V).\end{array}$$(17)

Emission-to-extinction ratios are subject to systematic errors because extinction probes the ISM in the foreground to the star, while emission probes the entire line of sight (Planck Collaboration Int. XXI 2015). Figure 22 presents the histogram of the ratio E(B − V)^⋆/E(B − V)^∞ of the reddening to the star to the total reddening, i.e., the fraction of ISM material that is in front of each star. The median ratio for our full sample is 0.83, illustrating the potential for systematic effects on the polarization ratios.

If we assumed for simplicity that the ISM along the line of sight were uniform (in density, magnetic-field orientation, and dust properties), then the polarization ratio R_P/p = P/p_V would artificially increase linearly with decreasing E(B − V)^⋆/E(B − V)^∞. Consequently, by neglecting the presence of background material we would typically overestimate the polarization ratio R_P/p by 17%. Given this contamination, to debias our estimate of R_P/p we replace p_V by a linear estimate of what its value would be if the star were at infinity:

$$\begin{array}{c}\hfill p_V^{\infty }=p_V\frac{E{(B-V)}^{\infty }}{E{(B-V)}^{\star }},\end{array}$$(18)

with an associated uncertainty

$$\begin{array}{c}\hfill {\sigma }_{p_V}^{\infty }={\sigma }_{p_V}\frac{E{(B-V)}^{\infty }}{E{(B-V)}^{\star }}.\end{array}$$(19)

We use similar expressions to estimate $q_V^{\infty }$ and $u_V^{\infty }$ and their uncertainties ${\sigma }_{u_V}^{\infty }$ and ${\sigma }_{q_V}^{\infty }$. On the other hand, R_S/V, as a ratio of non-dimensional quantities, would be unaffected by this uniform background.

6.3. Selection of the lines of sight

However, the ISM is not uniform and as a consequence our estimates of both R_S/V and R_P/p could be biased by the presence of a background whose properties are different from those of the foreground to the star. The magnitude of this uncertainty is evaluated in Appendix G.2. We minimize the contribution of this uncertainty by excluding those lines of sight with an important background, as inferred from the ratio E(B − V)^⋆/E(B − V)^∞ shown in Fig. 22. We explicitly choose to keep only stars for which E(B − V)^⋆/E(B − V)^∞ >  0.75.

We also exclude those lines of sight where Δψ_S/V is significantly different than the expected 0°, having found that for these lines of sight the rms scatter about the best fit correlations yielding R_P/p and R_S/V (Sect. 6.4) is indeed much larger. To be conservative and retain enough lines of sight for our subsequent statistical analysis, we excluded only lines of sight for which Δψ_S/V >  45°, i.e., about 9% of the sample. Our results are not sensitive to this particular choice.

Our final sample contains 1505 stars. The ISM towards these stars in emission is representative of diffuse dust at high Galactic latitudes, with MBB fit parameters T_d = 19.7 ± 1.3 K, β = 1.60 ± 0.15, and E(B − V)_τ ∈ [0.01, 0.24], with a mean reddening ⟨E(B − V)_τ⟩ = 0.03.

6.4. Determination of the polarization ratios

Following the approach in Planck Collaboration Int. XXI (2015), we derive R_P/p through a joint correlation of the pair (Q, U) with ($q_V^{\infty },u_V^{\infty }$), and derive R_S/V through a joint correlation of (Q/I, U/I) with (q_V/τ_V, u_V/τ_V). In Fig. 23 we present the two correlation scatterplots, that for R_P/p on the left, and that for R_S/V on the right. For the fitting, we settle on the Bayesian method of Kelly (2007), but we obtain the same results with other fitting methods making use of uncertainties on both axes. Both for determining the value of the ratio (the slope) and for calculating the reduced χ² to assess the quality of the fit, it is important to assess all sources of uncertainty, as addressed in this section and in more detail in Appendix G.

As a first test, we fit the data for the 206 translucent lines of sight from Planck Collaboration Int. XXI (2015) with this estimator and find no change, even though we smooth Planck Stokes parameter maps to 40′ FWHM. We are therefore confident that it is legitimate to compare the polarization ratios that we derive here at 40′ resolution with those measured at 7′ resolution in Planck Collaboration Int. XXI (2015).

The correlation of (Q, U) with ($q_V^{\infty },u_V^{\infty }$) shown in Fig. 23 (left) is tight, with a Pearson correlation coefficient −0.92. For lines of sight where p is low, error bars have been greatly increased by the correction factor $\sqrt{1+\delta p_{\mathrm{beam}}^2/p^2}$ for beam depolarization (Appendix G.1). When systematic uncertainties are taken into account, the reduced χ² is good (1.02, compared to 3.9 when they are ignored). The fit yields a polarization ratio R_P/p = (5.42 ± 0.05) MJy sr⁻¹, similar to the value found for translucent lines of sight, (5.4 ± 0.5) MJy sr⁻¹ (Planck Collaboration Int. XXI 2015).

We find a somewhat larger scatter in the correlation plot (Q/I, U/I) with (q_V/τ_V, u_V/τ_V) in Fig. 23 (right), as quantified by the slightly lower absolute value of the correlation coefficient. Here also, the reduced χ² is good when systematic uncertainties are included (0.96), and larger when they are not (3.0). The fitted slope R_S/V = 4.31 ± 0.04 is also compatible with the value 4.2 ± 0.5 found for translucent lines of sight at 7′ resolution (Planck Collaboration Int. XXI 2015). We find similar values at other resolutions: R_S/V = 4.2 at 80′ resolution; and R_S/V = 4.2 at 20′ resolution. These results are for the fiducial intensity offset. With the low and high offsets at 40′ resolution, we obtain R_S/V = 4.8 and R_S/V = 3.9, respectively, which makes the Planck intensity offset the main source of uncertainty on R_S/V.

6.5. Variations of R_P/p and R_S/V with column density

Our sample contains enough lines of sight to study variations of the polarization ratios with column density N_H, which is determined from the dust optical depth at 353 GHz as explained at the beginning of Sect. 5. This is presented in Fig. 24. The sample is divided into 15 independent bins in N_H, each bin containing 100 stars. For the polarization ratio R_P/p at low N_H, we observe a roughly 10% increase with increasing N_H, from about 5.0 MJy sr⁻¹ at N_H ≈ 6 × 10¹⁹ cm⁻² to 5.4 MJy sr⁻¹ at N_H ≈ 1.5 × 10²⁰ cm⁻². This is followed by a plateau at higher N_H. The normalized polarization ratio, R_S/V, is found to decrease with column density for the low total intensity offset, to be rather flat for the fiducial one, and to slightly increase for the high offset. The values obtained for R_P/p and R_S/V at higher (20′) and lower (80′) resolutions are close to the 40′ values. Both the Pearson correlation coefficients and the correlation coefficients provided by the Bayesian method (Kelly 2007) are high enough at all column densities to bring confidence in our results. For comparison, the values of the polarization ratios found for translucent lines of sight in Planck Collaboration Int. XXI (2015), together with their ranges of uncertainty, are also displayed in Fig. 24. Altogether, R_P/p and R_S/V are remarkably constant with column density and consistent with the values determined for translucent lines of sight. We note that there are some small variations of potential significance, such as the 10% increase in R_P/p over the range 6 × 10¹⁹ − 1.5 × 10²⁰ cm⁻², but it is beyond the scope of this paper to discuss them.

6.6. Maximum p_V/E(B − V) at low column densities

Regarding the observed starlight polarization fraction per unit reddening in the optical, p_V/E(B − V), its maximum value was first estimated to be at most 9% by Serkowski et al. (1975). This maximum value is often considered as providing an upper limit on the dust alignment efficiency, although it is based on less than 300 stars at moderate extinction (E(B − V) >  0.15), characteristic of translucent lines of sight. Several attempts have been made to constrain this maximum value at low reddening (Fosalba et al. 2002; Frisch et al. 2015; Skalidis et al. 2018), suggesting larger values, but with poor precision. Our present sample probes more diffuse lines of sight and, with more statistical significance, allows us to characterize the maximum polarization fraction at low reddening, extending to E(B − V) ≈ 0.01.

To study the dependence of the polarization fractions P/I and p_V/E(B − V) on column density, we combine data for the 1505 diffuse lines of sight to stars from this study (high latitude, low E(B − V), 40′ resolution for Planck data and E(B − V) map) and the 206 stars on translucent lines of sight (low latitude, moderate E(B − V), 7′ resolution for Planck data) from Planck Collaboration Int. XXI (2015).

Figure 25 shows how these polarization fractions in emission and extinction vary with column density. Polarization fractions at 353 GHz never reach low values because Planck polarization data have been corrected for beam depolarization through Eq. (F.63). However, this does not affect our results for the high percentiles because high values of p suffer very little depolarization. We overplot the upper limit p_V ≤ 9%E(B − V) of Serkowski et al. (1975), the non-linear fit p_V/E(B − V) ∝ E(B − V)^−0.2 proposed by Fosalba et al. (2002) for polarization in extinction (Fig. 25, right), and the estimate for the maximum value of the polarization fraction in emission observed for translucent lines of sight on the left (approximately 14%, Planck Collaboration Int. XXI 2015).

Assuming a maximum polarization in emission of p_max = 20% for our sample at a resolution of 40′ (close to its 99th percentile¹⁸), and a polarization ratio R_S/V = 4.3 (Fig. 23), we would expect a maximum polarization fraction in extinction of p_V/E(B − V)_max = (3.1 × p_max)/(R_S/V × 1.086)≈13%. As seen in Fig. 25 (right), this upper limit is somewhat smaller than the measured 99th percentile of our data in extinction. We would reach a similar conclusion using the value p_max = 22% that we obtained for the diffuse ISM in Sect. 4.1.1, with p_V/E(B − V)_max ≃ 14.5%. This upper limit is also smaller than the upper limit proposed by Fosalba et al. (2002) based on a study of the dependence of the mean starlight polarization fraction p_V/E(B − V) with E(B − V).

However, the distribution in the density of points for polarization in extinction (right panel) compared to that for polarization in emission (left panel) suggests that lines of sight with high p_V/E(B − V) might be outliers. One should indeed be aware that the direct derivation of the maximum polarization fraction in extinction is much more subject to noise and systematics than our derivation, which is based on the measurement of the polarization ratio R_S/V and the much better characterized maximum polarization fraction p in emission. We therefore consider the value of 13% as a well-constrained maximum value for p_V/E(B − V) at very low N_H (< 5 × 10²⁰ cm⁻²). This is a strong new constraint on the grain optical properties used in dust models.

The observed maximum polarization fractions drop from the diffuse ISM at high Galactic latitudes to the translucent lines of sight at low Galactic latitudes. In emission, p_max decreases from 20% to 14%, whereas in extinction (p_V/E(B − V))_max decreases from 13% to 9%. Such a decrease is usually attributed to a loss of grain alignment (see Andersson et al. 2015 and references therein). However, inspection of Fig. 21 for Galactic latitudes higher than 10° shows that the product 𝒮 × p remains constant over the range of column densities probed here. Following our analysis in Sect. 5, we therefore attribute most of this decrease in the maximum polarization fraction when the column density increases to the increasing depolarizing effect from the structure of the magnetic field along the line of sight, with little room for a systematic decrease in the grain alignment efficiency.

Dust models should therefore be able to reproduce the maximum observed polarization fractions, p_max = 20% in emission and (p_V/E(B − V))_max = 13% in extinction, even when applied to regimes in column densities where such values are actually never directly observed.

7. Conclusions

In this paper, we have analysed the Planck 2018 thermal dust polarization data at 353 GHz. Starting from full-sky maps of Stokes I, Q, and U at a uniform 80′ resolution, processed with the Generalized Needlet Internal Linear Combination (GNILC) algorithm (Remazeilles et al. 2011) to mitigate the contamination by CIB and CMB anisotropies as well as noise, we have presented the maps of polarization fraction p, polarization angle ψ, and polarization angle dispersion function 𝒮, with their associated uncertainties. The statistical analysis of these maps provides new insights into the astrophysics of dust polarization.

We have been able to determine the maximum observed polarization fraction, $p_{\text{max}}=22.0_{-1.4}^{+3.5}\pm 0.1\%$, at this resolution and frequency, where the second uncertainty is statistical, underscoring the excellent quality of the data, and the first reflects the principal systematic uncertainty affecting this determination, which is linked to the uncertainty on the Galactic emission zero level in total intensity (Planck Collaboration XI 2014). This maximum polarization fraction provides strong constraints for models of dust properties and alignment in the Galactic magnetic field (Guillet et al. 2018). Owing to the strong effect of the magnetic field morphology, a low value for the maximum polarization fraction in a given region is not an indication that grain alignment in that region is ineffective, but rather that polarization is strongly affected by depolarization because the direction of the large-scale field is closer to the line of sight in that region.

We confirmed that the statistical properties of p, ψ, and 𝒮 essentially reflect the structure of the Galactic magnetic field (Planck Collaboration Int. XIX 2015), with a clear peak of the polarization angle near 0°, corresponding to the field being parallel to the Galactic plane, and an inverse proportionality between the polarization fraction p and the polarization angle dispersion function 𝒮. We showed analytically, and using numerical models of the polarized sky, that this relationship can be reproduced statistically to first order by an interstellar magnetic field incorporating a turbulent component superposed on a small number of layers with a simple uniform field configuration. Looking for evidence in the diffuse ISM (N_H <  8 × 10²¹ cm⁻²) of a correlation of the polarization fraction with the dust temperature, as one could expect from the classical radiative torque theory (Lazarian & Hoang 2007), we could not find any: all variations of p are here again mirrored with those of 𝒮, which does not depend on the dust physics.

Based on this analysis, we introduced the product 𝒮 × p as a means of exploring the non-geometric elements of the polarization maps, such as variations in grain properties, in alignment physics, or in ISM phase contributions. We showed that 𝒮 × p exhibits smaller and smoother variations than either p or 𝒮 when considered as a function of the Galactic latitude b, the Galactic longitude l, or the column density (which is simply scaled from the dust optical depth τ at 353 GHz).

We provided an analysis at a finer angular resolution of 40′ using the Planck 2018 data, towards a sample of six molecular regions in the Gould Belt. This confirmed the trends observed at coarser resolution over the full sky, most notably that the polarization angle dispersion function is inversely proportional to the polarization fraction, 𝒮 ∝ 1/p. Strikingly, the 𝒮 versus p curves for the different regions all fall on the same line, demonstrating that most of the variations of p with column density are driven by variations of 𝒮, i.e., by the structure of the magnetic field along the line of sight (Planck Collaboration Int. XX 2015). Considering then the product 𝒮 × p and how it varies with dust temperature T_d in these regions, we found no evidence for a link between the polarization properties and the intensity of the radiation field. Comparing these properties with column density in a multi-resolution view, we found that the product 𝒮 × p does decrease, but only by about 25% between N_H ≈ 10²⁰ cm⁻² and N_H ≈ 2 × 10²² cm⁻², while the polarization fraction p decreases by a factor of 3−4 over the same interval.

We also compared the Planck polarization data with optical stellar polarization data in the high Galactic latitude sky, expanding on the analysis done in Planck Collaboration Int. XXI (2015) for translucent lines of sight. We constrained the polarization properties of the dust at low column densities, quantified by the polarization ratios R_P/p = P/p_V and R_S/V = (P/I)/(p/τ_V) defined in Planck Collaboration Int. XXI (2015). The larger statistical sample (1505 stars selected) allowed us to study the dependence of these polarization ratios on column density. We found R_P/p to increase from 5.0 MJy sr⁻¹ at N_H = 6 × 10¹⁹ cm⁻² to 5.4 MJy sr⁻¹ for N_H >  1.5 × 10²⁰ cm⁻², the same value as for translucent lines of sight. The polarization ratio R_S/V was found to be compatible on average (around 4.3) with that for translucent lines of sight (4.2 ± 0.5), having a decreasing, flat, or slightly increasing trend with the column density, depending on the offset for the Planck intensity map at 353 GHz, which is here again the dominant systematic effect. Combining emission and extinction measurements, we derived an estimate for the maximum value of the polarization fraction in the visible at high Galactic latitude, p_V/E(B − V) ≤ 13%, significantly higher than the value of 9% characterizing translucent lines of sight at low latitudes (Serkowski et al. 1975). This is a strong new constraint that dust models must satisfy.

---

Acknowledgments

The Planck Collaboration acknowledges the support of ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 Research & Innovation Framework Programme/ERC grant agreement ERC-2016-ADG-742719. This research has received funding from the Agence Nationale de la Recherche (ANR-17-CE31-0022). We thank Pekka Teerikorpi and Andrei Berdyugin for kindly providing stellar polarization data and providing insights on stellar polarization references, Gina Panopoulou for discussions of the zero point calibration of the polarization angle, and Ralf Siebenmorgen and Nikolai Voshchinnikov for statistical discussions. We gratefully acknowledge the help of Rosine Lallement regarding the handling of the Gaia data.

References

Appendix A: A guide to Planck papers on Galactic astrophysics using polarized thermal emission from dust

In this appendix, we give a summary of the contents and main results of the Planck papers dealing with Galactic astrophysics using polarized thermal emission from dust, in the hope that it will provide a useful reference for many readers.

In Planck Collaboration Int. XIX (2015), we presented the first analysis of the 353-GHz polarized sky at 1° resolution, focusing on the statistics of the polarization fraction p and polarization angle ψ, at low and intermediate Galactic latitudes. We found a high maximum polarization fraction (p_max >  19.8%) in the most diffuse regions probed. This maximum polarization fraction was found to decrease as total gas column density N_H increases, which might be interpreted as changes of grain alignment properties or as the effect of magnetic-field structure along the line of sight. We also characterized the structure of the polarization angle map by computing its local dispersion function 𝒮, which was found to be inversely related to the polarization fraction, lending support to the second explanation.

In Planck Collaboration Int. XX (2015), we presented an analysis of Planck 353-GHz polarized thermal dust emission towards a set of molecular clouds and other nearby fields, at 15′ resolution, and compared their statistics to those computed on synthetic maps derived from simulations of anisotropic magnetohydrodynamical (MHD) turbulence. We showed that, at these angular scales, the turbulent structure of the Galactic magnetic field (GMF) is able to reproduce the main statistical properties of polarized thermal dust emission in nearby molecular clouds, with no necessity to introduce spatial changes of dust alignment properties.

In Planck Collaboration Int. XXI (2015), we compared the Planck polarized emission at 353 GHz to surveys of starlight polarization in extinction in the visible, selecting those stars for which the submillimetre and optical estimates of column densities and polarization angles match. For these lines of sight, we computed the ratio R_S/V of submillimetre to visible polarization fractions, and the ratio R_P/p of the polarized intensity in the submillimetre to the degree of polarization in the visible. We found these to be R_S/V = 4.2 ± 0.2 ± 0.3 and R_P/p = [5.4 ± 0.2 ± 0.3] MJy sr⁻¹, where the first uncertainty is statistical and the second is systematic. The value of R_P/p provides strong constraints for models of dust polarized emission. The DustEM model (Compiègne et al. 2011) has been updated by Guillet et al. (2018) to take these constraints into account.

In Planck Collaboration Int. XXXII (2016), we studied the correlation between the magnetic-field orientation inferred from polarization angles at 353 GHz and the filamentary structures of matter, at 15′ resolution, for intermediate and high Galactic latitudes, covering a range of total gas column densities from 10²⁰ cm⁻² to 10²² cm⁻². The filaments were extracted using the Hessian matrix of the dust total intensity map. We found that the filaments are preferentially parallel to the magnetic orientation, in particular when the polarization fractions are high and the filaments are more diffuse. Conversely, some of the densest, molecular filaments are perpendicular to the magnetic orientation. This analysis also provided a first estimate for the ratio of the turbulent to mean components of the GMF, i.e., f_M = 0.8 ± 0.2.

In Planck Collaboration Int. XXXIII (2016), we further studied the signature of the magnetic-field geometry of interstellar filaments in Planck 353-GHz dust polarization maps, at the native 4.​​′8 resolution, focusing on three nearby, dense, star-forming filaments (N_H ≈ 10²² cm⁻²), and interpreting the Planck Stokes I, Q, and U maps as the superposition of contributions from the filaments themselves and their respective backgrounds and foregrounds. In this way we found differences in polarization angles between the filaments and their environments, reaching values up to 54°, and a decrease in polarization fraction within the filaments, although this could be due not only to the effect of magnetic field tangling along the line of sight, but also in part to changes of grain alignment properties.

In Planck Collaboration Int. XXXIV (2016), we combined the polarization data from Planck at 353 GHz with rotation measure (RM) observations from Savage et al. (2013) towards a massive star-forming region, the Rosette Nebula in the Monoceros molecular cloud, to study the impact of an expanding H II region on the morphology of the Galactic magnetic field. We introduced an analytical model that describes the magnetic field structure in a spherical shell, following the expansion of an ionized nebula in a medium with uniform density and magnetic field, and fitted it to the data. This work was subsequently extended to non-spherical bubbles to model the structure of the large-scale magnetic field in the Local Bubble (Alves et al. 2018). The Planck polarization data towards the Orion–Eridanus superbubble provide additional evidence for the impact of massive stars on the magnetic field structure (Soler et al. 2018).

In Planck Collaboration Int. XXXV (2016), we studied the relative orientation between filamentary structures of matter and the magnetic field towards molecular clouds of the Gould Belt, probing lines of sight with total gas column densities N_H from around 10²¹ to 10²³ cm⁻², at 10′ resolution, using the histogram of relative orientations (HRO) technique (Soler et al. 2013). We found that this relative orientation changes progressively from preferentially parallel in regions with the lowest gas column densities to preferentially perpendicular in the regions with the highest N_H, with a crossover at N_H of a few 10²¹ cm⁻². This change in relative orientation was found to be compatible with simulations of sub- or trans-Alfvénic MHD turbulence, underlining the important dynamical role played by the magnetic field in shaping the structure of molecular clouds.

In Planck Collaboration Int. XXXVIII (2016), we studied the E and B modes (Kamionkowski et al. 1997; Zaldarriaga & Seljak 1997) of dust polarization from the magnetized filamentary structure of the interstellar medium at high Galactic latitudes, isolating Hessian-extracted filaments at angular scales where the E/B power-asymmetry and TE correlation are observed (Planck Collaboration Int. XXX 2016). The preferred orientation of these filaments parallel to the magnetic field is able to account for both of these observations and was quantified by an oriented stacking of the maps of I, Q, U, E, and B. From these stacked maps and the histogram of relative orientations, we derived an estimate of the mean polarization fraction in the filaments, ⟨p⟩≈11%.

In Planck Collaboration Int. XLII (2016), we provided a comparison of three different models of the large-scale GMF to Planck polarization data at low and high frequencies, respectively taken as templates for the polarized synchrotron and thermal dust emission. We found in particular that the models underpredict the dust polarization out of the Galactic plane, calling for an increased ordering of the GMF near the observer.

In Planck Collaboration Int. XLIV (2016), we provided a phenomenological model of the polarized dust sky. Polarized emission is assumed to arise from a small number of layers in which the GMF is taken to be a superposition of a uniform component B₀ and a turbulent component B_t. Applying this model to the Planck maps of the southern Galactic cap at 353 GHz and 1° resolution, and using the 1-point statistics of p and ψ, we could constrain the orientation of the large-scale field in the Solar neighbourhood, the number of layers (N ≈ 7), the effective polarization fraction of dust emission (p₀ = 26 ± 3%), and the ratio of the strengths of the turbulent to mean components of the GMF (f_M = 0.9 ± 0.1). This phenomenological framework was further improved by Ghosh et al. (2017) and Vansyngel et al. (2017).

Appendix B: GNILC Stokes and covariance maps

B.1. Variable resolution GNILC maps

For reference, in Fig. B.1 we show the GNILC Stokes maps at 353 GHz and variable resolution over the sky, alongside the map of the effective FWHM, whose discrete values are 5′, 7′, 10′, 15′, 20′, 30′, 60′, and 80′. The total intensity map corresponds to the fiducial offset value.

B.2. GNILC-processed covariance maps

To assess the statistical uncertainties affecting the dust polarization properties studied in this paper, the GNILC algorithm provides the maps of covariances between the Stokes parameters at 353 GHz, σ_II, σ_IQ, σ_IU, σ_QQ, σ_QU, and σ_UU. We now describe how these covariance maps are produced.

The GNILC dust map, D_GNILC, is a mimimum-variance weighted linear combination of the Planck frequency maps X_i:

$$\begin{array}{c}\hfill D_{\mathtt{GNILC}}={\displaystyle \sum _i{w}_i{X}_i,}\end{array}$$(B.1)

where the sum extends over the seven Planck polarization-sensitive frequency channels. The weights w_i are estimated by the GNILC algorithm in order to extract the dust emission while filtering out the instrumental noise and the CMB¹⁹. The residual noise rms fluctuations, N_GNILC, of the GNILC dust map are thus related to the instrumental noise rms fluctuations, N_i, in each Planck frequency map X_i, according to

$$\begin{array}{c}\hfill N_{\mathtt{GNILC}}={\displaystyle \sum _i{w}_i{N}_i.}\end{array}$$(B.2)

That residual noise is minimized by the GNILC weights. An estimate of the instrumental noise rms fluctuations N_i in each frequency channel i can be obtained by computing the half-difference of the Planck half-mission maps X_i, HM1 and X_i, HM2:

$$\begin{array}{c}\hfill {\widehat{N}}_i=\frac{1}{2}(X_{i,{\mathrm{HM}}_2}-X_{i,{\mathrm{HM}}_1}),\end{array}$$(B.3)

because the sky emission cancels out in the difference while the noise does not. We can thus estimate the residual noise in the GNILC dust map as

$$\begin{array}{c}\hfill {\widehat{N}}_{\mathtt{GNILC}}={\displaystyle \sum _i{w}_i{\widehat{N}}_i,}\end{array}$$(B.4)

where the ${\widehat{N}}_i$ maps have first been smoothed to the actual resolution of the GNILC maps, i.e., either 80′ for the uniform resolution case or to the local resolution of the specific regions of the sky for the multi-resolution case. The estimate ${\widehat{N}}_{\mathtt{GNILC}}$ has the variance of the actual residual noise in the GNILC dust maps, D_GNILC.

The GNILC noise covariance maps are then estimated as follows. We first smooth the native Planck noise covariance maps at 353 GHz, σ_XY²⁰, where X and Y stand for any one of the three Stokes I, Q, or U, to the resolution of the GNILC maps, by following the procedure employed in Planck Collaboration Int. XIX (2015). For the multi-resolution case, the value of the smoothing scale adopted in each region of the sky depends on the local resolution of the GNILC maps in that region. Because a covariance is derived from the product of two Stokes parameter maps, the smoothing scale of a covariance map is $\sqrt{2}$ times the resolution of the Stokes maps. We then compute the local (co)variance value in each region of a given resolution of the Planck and GNILC noise maps, ${\widehat{N}}_{353}$ and ${\widehat{N}}_{\mathtt{GNILC}}$, for instance:

$${\text{cov}}^{(j)}{\widehat{N}}_{\text{353}}^Q,{\widehat{N}}_{\text{353}}^U=\frac{1}{n_j}{\displaystyle \sum _{p\in R_j}}{\widehat{N}}_{\text{353}}^Q(p){\widehat{N}}_{\text{353}}^U(p),$$(B.5)

$$\begin{array}{cc}\hfill & {\mathrm{cov}}^{(j)}({\widehat{N}}_{\mathtt{GNILC}}^Q,{\widehat{N}}_{\mathtt{GNILC}}^U)=\frac{1}{n_j}{\displaystyle \sum _{p\in R_j}{\widehat{N}}_{\mathtt{GNILC}}^Q(p){\widehat{N}}_{\mathtt{GNILC}}^U(p),}\hfill \end{array}$$(B.6)

where R_j is the set of sky pixels at the given GNILC effective resolution, indexed by j, and n_j is the number of such pixels. In each region R_j, the Planck noise covariance maps, σ_XY(p), are then scaled according to

$$\begin{array}{c}\hfill {\sigma }_{XY,\mathtt{GNILC}}(p)=\frac{{\mathrm{cov}}^{(j)}({\widehat{N}}_{\mathtt{GNILC}}^X,{\widehat{N}}_{\mathtt{GNILC}}^Y)}{{\mathrm{cov}}^{(j)}({\widehat{N}}_{353}^X,{\widehat{N}}_{353}^Y)}{\sigma }_{\mathit{XY}}(p).\end{array}$$(B.7)

The resulting covariance maps σ_{XY, GNILC} (p) are what we refer to as the GNILC-processed covariance maps in the rest of this paper.

We note that they are built using the PSB-only data for polarization at 353 GHz, both for the uniform 80′ resolution case and for the B-mode-driven, varying resolution case (5′–80′). In both cases, they are computed at a HEALPix resolution N_side = 2048, but in the uniform 80′ resolution case, they are downgraded to N_side = 128 to avoid oversampling. For the varying resolution case, we keep the original N_side = 2048. The maps are also converted from ${\mathrm{K}}_{\mathrm{CMB}}^2$ to MJy² sr⁻² using the conversion factor $287.5\mathrm{MJy}{\mathrm{sr}}^{-1}{\mathrm{K}}_{\mathrm{CMB}}^{-1}$ at 353 GHz (Planck Collaboration III 2020).

Figure B.2 shows these covariance maps for the variable resolution case, while Fig. B.3 shows the covariance maps at the common, uniform 80′ resolution. The sky patterns of these uniform resolution covariance maps are by construction extremely similar to those directly taken from the Planck 2018 data, but with a significant improvement in the amplitudes, as shown in Table B.1, which gives the mean ratios σ_XY, GNILC/σ_XY of the GNILC covariance maps σ_XY,GNILC at 80′ resolution and N_side = 128 to the corresponding maps σ_XY in the Planck 2018 data release. Consequently, we use these GNILC-processed covariance maps to assess statistical uncertainties in our analysis.

Appendix C: End-to-end simulations

The quality of the data presented here is assessed through a series of end-to-end (E2E) simulations that take into account all the known systematics and noise properties of the data. The process begins with a model of the sky (including foregrounds, in both total intensity and polarization), from which timelines (including all known effects) are simulated. These timelines are then processed through the Planck mapmaking pipeline (Planck Collaboration III 2020). These E2E simulations are the same ones used in Planck Collaboration XI (2020) and we refer the reader to appendix A of that paper for a detailed description. The dust component of the model used in these simulations is a combination of a realization of the Vansyngel et al. (2017) statistical model for ℓ ≥ 20 and the actual Planck 2018 Q and U maps at 353 GHz for ℓ ≤ 10, with a smooth transition between the two in the 10 ≤ ℓ ≤ 20 range. Because of the latter, the large-scale component of the field varies over the simulated sky, which is essential to reproduce the statistics of p. The other input sky components are taken from the latest version of the Planck Sky Model (Planck Collaboration XII 2016). These model component maps are then combined with the first 100 realizations of the systematic effects and noise (Planck Collaboration III 2020). This results in 100 E2E I, Q, and U maps at 353 GHz, which we smooth to 60′ resolution, and from which we derive polarization quantities and compare them to the input dust maps, after subtraction of the CMB and the CIB monopole, as is done for the ASMs (Sect. 2.4). This allows us to assess the effects of residual systematics and data noise on the statistics presented in this paper.

Figure C.1 shows that the difference between the input dust polarization fraction p⁽⁰⁾ and the average ⟨p⟩ over the 100 realizations of the E2E simulations is at most around 1%. The distribution function of the difference p − p⁽⁰⁾ is peaked around zero for each simulation. The average of these distributions is shown in the bottom panel of Fig. C.1, along with the 1σ dispersion around the average distribution, which has a mean of 0.03% and a standard deviation of 0.47%.

For the polarization angle data, the diagnostic of the E2E simulations is shown in Fig. C.2. To compute the average difference between the output and input polarization angles, ⟨ψ⟩−ψ⁽⁰⁾, account is taken of the circularity of the difference for each sky pixel and each realization independently. One can see that the regions of the sky where this difference is the largest are those where ψ⁽⁰⁾ crosses the ±90° boundary. The average distribution function of these differences over the 100 simulations and the 1σ dispersion about the average are shown in the bottom panel of the figure. The average distribution has a mean of 0${}_{.}^{\circ }$3 and a standard deviation of 8${}_{.}^{\circ }$3.

The same diagnostics are run on the polarization angle dispersion function 𝒮, which at 60′ resolution we compute with a lag δ = 30′. Results are shown in Fig. C.3. We note that the average ⟨𝒮⟩ of the output 𝒮 maps exhibits a significant positive bias with respect to the input 𝒮⁽⁰⁾ map, especially towards the Galactic poles. The distributions of pixel values in difference maps 𝒮 − 𝒮⁽⁰⁾ from the 100 simulations have a positive skewness. This shows that the polarization angle dispersion function is still affected by residual bias at this resolution, even though it is barely affecting the polarization angle map itself.

For completeness, Fig. C.4 shows the histograms of the polarization fractions, polarization angles, and polarization angle dispersions for the input (black curve), and the outputs of the E2E simulations. For the latter, the blue curve on each panel shows the average histogram over the 100 simulations, with the ±1σ dispersion among histograms shown as the blue area between dashed lines. The agreement is excellent for both quantities p and ψ, but for 𝒮 we note that at intermediate values the positive bias already mentioned appears clearly. Finally, we stress that although the polarization fractions rarely go above 20% for these simulated dust maps, this does not mean that the same range is expected in the Planck data.

Appendix D: Link between 𝒮 and the polarization gradients

D.1. Analytical derivation

The link between the polarization angle dispersion function 𝒮 and the polarization angle gradient |∇ψ| can be established analytically via a Taylor expansion of the polarization angle difference appearing in the definition of 𝒮:

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\mathit{r},\delta )=\langle {[\psi (\mathit{r}+\mathit{\delta })-\psi \left(\mathit{r}\right)]}^2\rangle \approx \langle {[\mathit{\delta }.\mathbf{\nabla }\psi ]}^2\rangle ,\end{array}$$(D.1)

where the average is computed over the annulus centred on r having inner and outer radii δ/2 and 3δ/2, respectively (Eq. (4)), and ∇ψ is the vector gradient of the polarization angle at the centre r. Using a local reference frame with axes y and z in the plane of the sky, we can write the displacement vector as

$$\begin{array}{c}\hfill \mathit{\delta }=l(cos\theta {\mathit{e}}_y+sin\theta {\mathit{e}}_z),\end{array}$$(D.2)

with δ/2 ≤ l ≤ 3δ/2. The expression of 𝒮² therefore becomes

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\mathit{r},\delta )\approx \langle l^2\rangle \langle {(cos\theta \frac{\partial \psi }{\partial y}+sin\theta \frac{\partial \psi }{\partial z})}^2\rangle ,\end{array}$$(D.3)

where the spatial average takes into account that l and θ are independent variables. The former simply yields

$$\begin{array}{c}\hfill \langle l^2\rangle =\frac{13}{12}{\delta }^2\approx {\delta }^2\end{array}$$(D.4)

and the latter average is over θ ∈ [0, 2π]. In that average, taking the square and averaging over θ cancels the cross product, so that

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\mathit{r},\delta )\approx \frac{{\delta }^2}{2}[{\left(\frac{\partial \psi }{\partial y}\right)}^2+{\left(\frac{\partial \psi }{\partial z}\right)}^2].\end{array}$$(D.5)

On the other hand, by defining the angular polarization gradient (Eq. (6)) for a unit polarization vector Q/P = cos(2ψ) and U/P = sin(2ψ), we have

$$\begin{array}{c}\hfill |\mathrm{\nabla }\psi |=2\sqrt{{\left(\frac{\partial \psi }{\partial y}\right)}^2+{\left(\frac{\partial \psi }{\partial z}\right)}^2},\end{array}$$(D.6)

which leads to the relation

$$\begin{array}{c}\hfill \mathcal{S}(\mathit{r},\delta )\approx \frac{\delta }{2\sqrt{2}}|\mathrm{\nabla }\psi |.\end{array}$$(D.7)

D.2. The case of Planck data

Figure D.1 shows the maps of both polarization gradients, |∇ψ| and |∇P| from Eqs. (6) and (5), respectively, for the GNILC-processed Planck data at 353 GHz and 160′ resolution. The correlations between |∇ψ| and 𝒮 on the one hand, and between |∇P| and 𝒮 on the other, are shown in Fig. D.2 (for 𝒮 a lag of 80′ is used). These plots show that 𝒮 correlates well with |∇ψ|, but not as well with |∇P| and that |∇ψ| is a very good proxy for the angular dispersion function 𝒮, as would be expected from Eq. (D.7) (and much faster to compute in practice).

Appendix E: Polarization fraction versus total gas column density for low and high offsets

In this appendix, we show in Fig. E.1 plots similar to Fig. 9, but for the low and high total intensity offsets (Sect. 2.2). The effects of the offset are discussed in Sect. 4.2.1.

Appendix F: The inverse relationship between 𝒮 and p

In this appendix, we use a phenomenological model of the submillimetre polarized thermal dust emission, developed in Planck Collaboration Int. XLIV (2016), Ghosh et al. (2017), and Vansyngel et al. (2017), to derive the relationship between the polarization fraction p and the polarization angle dispersion function 𝒮. In its most basic form presented in Fig. F.1, this model assumes the polarized sky to result from a small set of N concentric layers, each emitting a fraction 1/N of the total intensity²¹, and harbouring a magnetic field B = B₀ + B_t, where B₀ is a uniform field (the same in each layer) and B_t is an isotropic turbulent field that is taken, in each layer, as a different realization of a Gaussian random field in three dimensions. No effects of dust evolution or changes of intrinsic polarization properties of the dust grains are included in the model. By design, this model is able to reproduce the 1-point statistics of polarized thermal Galactic dust emission maps observed by Planck, but it turns out that it is also able to reproduce the trend 𝒮 ∝ 1/p and the probability density function of 𝒮 × p, as we demonstrate below.

F.1. Reference frame and notations

We use a reference frame defined in Fig. F.2. The x axis is the line of sight, oriented towards the observer, and Oyz is the plane of the sky. In that frame, the components of the large-scale magnetic field are (B_0x, B_0y, B_0z).

F.2. Magnetic field in a layer at a given line of sight

We begin by noting that, from one layer to the next, the different Gaussian realizations of the turbulent magnetic field can be taken to be independent. Therefore, in each layer i (with 1 ≤ i ≤ N), we can write the components of the magnetic field B_i = B₀ + B_t, at the position considered as the central pixel in the definition of 𝒮 (Eq. (4)), as

$$\begin{array}{cc}\hfill B_{\mathit{ix}}=& B_{0x}+f_{\mathrm{M}}{B}_0{G}_x,\hfill \end{array}$$(F.1)

$$\begin{array}{cc}\hfill B_{\mathit{iy}}=& B_{0y}+f_{\mathrm{M}}{B}_0{G}_y,\hfill \end{array}$$(F.2)

$$\begin{array}{cc}\hfill B_{\mathit{iz}}=& B_{0z}+f_{\mathrm{M}}{B}_0{G}_z.\hfill \end{array}$$(F.3)

Here G_x, G_y, and G_z are Gaussian random variables with zero mean and variance 1/3. The parameter

$$\begin{array}{c}\hfill f_{\mathrm{M}}=\frac{{\sigma }_{\mathrm{B}}}{B_0},\end{array}$$(F.4)

is the ratio of the standard deviation ${\sigma }_B=\sqrt{⟨{{\mathit{B}}_{\text{t}}}^2⟩}$ of the turbulent magnetic field to the magnitude B₀ = ||B₀|| of the ordered field. The orientation of the magnetic field at the central pixel in layer i is given by a set of two angles Γ_i ∈ [ − π/2, π/2] and ϕ_i ∈ [0, 2π]:

$$\begin{array}{cc}\hfill B_{\mathit{ix}}& =B_{i}sin{\mathrm{\Gamma }}_i;\hfill \end{array}$$(F.5)

$$\begin{array}{cc}\hfill B_{\mathit{iy}}& =B_{i}cos{\mathrm{\Gamma }}_{i}sin{\varphi }_i;\hfill \end{array}$$(F.6)

$$\begin{array}{cc}\hfill B_{\mathit{iz}}& =B_{i}cos{\mathrm{\Gamma }}_{i}cos{\varphi }_i.\hfill \end{array}$$(F.7)

As presented in Fig. F.2, the angle Γ_i is the inclination angle of the vector B_i with respect to the plane of the sky, while ϕ_i is the angle, counted positively clockwise from the north, between the z axis and the projection B_POS of B_i onto the plane of the sky.

F.3. Fluctuations within each layer over the scale δ

When computing the polarization angle dispersion function 𝒮, we introduce a specific scale, the lag δ, which we always take as half the FWHM ω, so that δ = ω/2. Presumably, the orientation of the magnetic field in each layer, i.e., the angles Γ_i and ϕ_i, vary little over these scales. Let us therefore consider a small Gaussian fluctuation δB_i around the direction of B_i. The rms σ_Bi(δ) of this fluctuation can be cast into a parameter similar in form to Eq. (F.4),

$$\begin{array}{c}\hfill f_{\mathrm{m}}(\delta )=\frac{{\sigma }_{B_i}(\delta )}{B_0}\ll 1,\end{array}$$(F.8)

which depends on the lag δ considered, and is related to the overall turbulence parameter f_M and to the spectral index α_M of the magnetic field²². This fluctuation δB_i corresponds to small variations for angles δΓ_i and δϕ_i:

$$\begin{array}{cc}\hfill \frac{\delta B_{\mathit{ix}}}{B_0}=& cos{\mathrm{\Gamma }}_i\delta {\mathrm{\Gamma }}_i=f_{\mathrm{m}}(\delta )g_x;\hfill \end{array}$$(F.9)

$$\begin{array}{cc}\hfill \frac{\delta B_{\mathit{iy}}}{B_0}=& -sin{\mathrm{\Gamma }}_{i}sin{\varphi }_i\delta {\mathrm{\Gamma }}_i+cos{\mathrm{\Gamma }}_{i}cos{\varphi }_i\delta {\varphi }_i=f_{\mathrm{m}}(\delta )g_y;\hfill \end{array}$$(F.10)

$$\begin{array}{cc}\hfill \frac{\delta B_{\mathit{iz}}}{B_0}=& -sin{\mathrm{\Gamma }}_{i}cos{\varphi }_i\delta {\mathrm{\Gamma }}_i-cos{\mathrm{\Gamma }}_{i}sin{\varphi }_i\delta {\varphi }_i=f_{\mathrm{m}}(\delta )g_z.\hfill \end{array}$$(F.11)

where g_x, g_y, and g_z are Gaussian random variables with zero mean and variance 1/3. This allows us to compute the small variations of the angles as

$$\begin{array}{cc}\hfill \delta {\mathrm{\Gamma }}_i& =\frac{f_{\mathrm{m}}(\delta )}{cos{\mathrm{\Gamma }}_i}{g}_x,\hfill \end{array}$$(F.12)

$$\begin{array}{cc}\hfill \delta {\varphi }_i& =\frac{f_{\mathrm{m}}(\delta )}{cos{\mathrm{\Gamma }}_i}(g_{y}cos{\varphi }_i-g_{z}sin{\varphi }_i),\hfill \end{array}$$(F.13)

provided that the ratio $\stackrel{\text{m}}{.}/cos{\mathrm{\Gamma }}_i$ is still small, which only fails if the direction of the mean magnetic field is very close to the line of sight.

F.4. Polarization angle and Stokes parameters

The angle ϕ_i is related to the polarization angle ψ_i, appearing in the definition of the Stokes parameters below, by a π/2 rotation, i.e., ψ_i = ϕ_i − π/2[π]. The π-ambiguity arises because the Stokes parameters are unchanged in the transformation B_POS ↦ −B_POS. The polarization angle thus lies in the range [ − π/2, π/2], and the Stokes parameters (I_i, Q_i, U_i) at the central pixel for each layer i are then²³

$$\begin{array}{cc}\hfill Q_i& =p_{\max }{I}_i{cos}^2{\mathrm{\Gamma }}_{i}cos2{\psi }_i=P_{i}cos2{\psi }_i,\hfill \end{array}$$(F.14)

$$\begin{array}{cc}\hfill U_i& =p_{\max }{I}_i{cos}^2{\mathrm{\Gamma }}_{i}sin2{\psi }_i=P_{i}sin2{\psi }_i,\hfill \end{array}$$(F.15)

where I_i and P_i are the total and polarized intensity at the central pixel in layer i, respectively, and p_max is the polarization fraction of thermal dust emission that would be observed in the case of a uniform magnetic field parallel to the plane of the sky (Γ_i = 0). A fluctuation δB_i of the magnetic field at the scale δ therefore produces a small variation of these Stokes parameters that can be written as

$$\begin{array}{cc}\hfill \delta Q_i=& -2(Q_{i}tan{\mathrm{\Gamma }}_i\delta {\mathrm{\Gamma }}_i+U_i\delta {\psi }_i),\hfill \end{array}$$(F.16)

$$\begin{array}{cc}\hfill \delta U_i=& -2(U_{i}tan{\mathrm{\Gamma }}_i\delta {\mathrm{\Gamma }}_i-Q_i\delta {\psi }_i),\hfill \end{array}$$(F.17)

where it is assumed that the total intensity varies little on the scale δ. Because we work with a lag smaller than the FWHM, δ = ω/2, this is a reasonable assumption. The fluctuation of the polarization angle is δψ_i = δϕ_i, and so by inserting Eqs. (F.12) and (F.13) for the fluctuations of the angles we obtain

$$\begin{array}{cc}\hfill \delta Q_i=& -\frac{2f_{\mathrm{m}}(\delta )}{cos{\mathrm{\Gamma }}_i}[Q_i{g}_{x}tan{\mathrm{\Gamma }}_i+U_i(g_{y}cos{\varphi }_i-g_{z}sin{\varphi }_i)],\hfill \end{array}$$(F.18)

$$\begin{array}{cc}\hfill \delta U_i=& -\frac{2f_{\mathrm{m}}(\delta )}{cos{\mathrm{\Gamma }}_i}[U_i{g}_{x}tan{\mathrm{\Gamma }}_i-Q_i(g_{y}cos{\varphi }_i-g_{z}sin{\varphi }_i)].\hfill \end{array}$$(F.19)

These expressions will be helpful in determining the fluctuations of the Stokes parameters over which to average when computing the polarization angle dispersion function in the next section.

F.5. Polarization angle dispersion function

The polarization angle dispersion function 𝒮 is computed for a central pixel c, and consists of an average over the n pixels, indexed by j (with 1 ≤ j ≤ n), in an annulus of mean radius δ = ||δ|| and width δ around the central pixel, as defined in Eq. (4). This can also be written in terms of the Stokes parameters Q and U at the central pixel, and Q(j) and U(j) at a pixel j in the annulus (Planck Collaboration Int. XIX 2015):

$$\begin{array}{c}\hfill \mathcal{S}(\delta )={\displaystyle \sqrt{\frac{1}{n}\sum _{j=1}^n{[\frac{1}{2}arctan\frac{Q(j)U-U(j)Q}{Q(j)Q+U(j)U}]}^2}.}\end{array}$$(F.20)

Because we are interested in the average behaviour of 𝒮, we will ultimately consider the mean of this expression over the position of the central pixel as well.

The distribution function of 𝒮 (Fig. 7) shows that most pixels have a small dispersion of polarization angles, 𝒮 ≲ 10°. For these values, it is safe to approximate the arctangent by its argument, so that we may write

$$\begin{array}{c}\hfill 4{\mathcal{S}}^2(\delta )={\langle {\left[\frac{Q(j)U-U(j)Q}{Q(j)Q+U(j)U}\right]}^2\rangle }_j.\end{array}$$(F.21)

The Stokes parameters at pixels c and j can be written as sums over the N layers. More precisely, for the central pixel we have, by definition,

$$\begin{array}{cc}\hfill Q=& {\displaystyle \sum _{i=1}^N{Q}_i,}\hfill \end{array}$$(F.22)

$$\begin{array}{cc}\hfill U=& {\displaystyle \sum _{i=1}^N{U}_i,}\hfill \end{array}$$(F.23)

while for the displaced pixel j we can write

$$\begin{array}{cc}\hfill Q(j)=& {\displaystyle \sum _{i=1}^N[Q_i+\delta Q_i(j)]=Q+\delta Q(j),}\hfill \end{array}$$(F.24)

$$\begin{array}{cc}\hfill U(j)=& {\displaystyle \sum _{i=1}^N[U_i+\delta U_i(j)]=U+\delta U(j),}\hfill \end{array}$$(F.25)

exhibiting the fluctuations of the Stokes parameters given in Eqs. (F.18) and (F.19). We use this decomposition to write the numerator and denominator that appear in the squared quantity above as

$$\begin{array}{cc}\hfill Q(j)U-U(j)Q=& U\delta Q(j)-Q\delta U(j),\hfill \end{array}$$(F.26)

$$\begin{array}{cc}\hfill Q(j)Q+U(j)U=& P^2+Q\delta Q(j)+U\delta U(j).\hfill \end{array}$$(F.27)

In the latter expression, the second and third terms are most likely negligible compared to the polarized emission P² at the central pixel, especially when averaged over index j, and so they can be ignored. We therefore have

$$\begin{array}{c}\hfill 4{\mathcal{S}}^2(\delta )\approx \frac{{\langle {[U\delta Q(j)-Q\delta U(j)]}^2\rangle }_j}{P^4},\end{array}$$(F.28)

because P is independent of the pixel j. Appearing in the numerator are the averages ⟨δQ(j)²⟩_j, ⟨δU(j)²⟩_j, and ⟨δQ(j)δU(j)⟩_j, which can be computed using the fact that for the above Gaussian random variables g_x, g_y, and g_z

$$\begin{array}{cc}& {\langle g_x^2\rangle }_j={\langle g_y^2\rangle }_j={\langle g_z^2\rangle }_j=\frac{1}{3},\hfill \end{array}$$(F.29)

$$\begin{array}{cc}& {\langle g_x{g}_y\rangle }_j={\langle g_x{g}_z\rangle }_j={\langle g_y{g}_z\rangle }_j=0.\hfill \end{array}$$(F.30)

Because the random variables g_x, g_y, and g_z are also uncorrelated from one layer to the next, we have

$$\begin{array}{cc}& {\langle \delta Q{(j)}^2\rangle }_j={\displaystyle \sum _{i=1}^N{\langle \delta Q_i{(j)}^2\rangle }_j,}\hfill \end{array}$$(F.31)

$$\begin{array}{cc}& {\langle \delta U{(j)}^2\rangle }_j={\displaystyle \sum _{i=1}^N{\langle \delta U_i{(j)}^2\rangle }_j,}\hfill \end{array}$$(F.32)

$$\begin{array}{cc}& {\langle \delta Q(j)\delta U(j)\rangle }_j={\displaystyle \sum _{i=1}^N{\langle \delta Q_i(j)\delta U_i(j)\rangle }_j,}\hfill \end{array}$$(F.33)

which yields, using the expressions of Eqs. (F.18) and (F.19),

$$\begin{array}{c}\hfill {\langle \delta Q{(j)}^2\rangle }_j=\frac{4}{3}{f}_{\mathrm{m}}^2(\delta ){\displaystyle \sum _{i=1}^N\frac{Q_i^2{tan}^2{\mathrm{\Gamma }}_i+U_i^2}{{cos}^2{\mathrm{\Gamma }}_i},}\end{array}$$(F.34)

$$\begin{array}{c}\hfill {\langle \delta U{(j)}^2\rangle }_j=\frac{4}{3}{f}_{\mathrm{m}}^2(\delta ){\displaystyle \sum _{i=1}^N\frac{U_i^2{tan}^2{\mathrm{\Gamma }}_i+Q_i^2}{{cos}^2{\mathrm{\Gamma }}_i},}\end{array}$$(F.35)

$$\begin{array}{c}\hfill {\langle \delta Q(j)\delta U(j)\rangle }_j=\frac{4}{3}{f}_{\mathrm{m}}^2(\delta ){\displaystyle \sum _{i=1}^N\frac{Q_i{U}_i({tan}^2{\mathrm{\Gamma }}_i-1)}{{cos}^2{\mathrm{\Gamma }}_i}.}\end{array}$$(F.36)

Combining the above expressions, we then have

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\delta )=\frac{f_{\mathrm{m}}^2(\delta )}{3P^4}{\displaystyle \sum _{i=1}^N\frac{{(Q{Q}_i+U{U}_i)}^2+{(Q{U}_i-U{Q}_i)}^2{tan}^2{\mathrm{\Gamma }}_i}{{cos}^2{\mathrm{\Gamma }}_i}.}\end{array}$$(F.37)

The combinations of Q, U, Q_i, and U_i appearing in this expression can be cast into another form by introducing the angular shift Δψ_i = ψ_i − ψ between the polarization angle in each layer ψ_i and the observed polarization angle ψ, both considered at the central pixel:

$$\begin{array}{cc}\hfill Q{Q}_i+U{U}_i=& P_{i}Pcos2\mathrm{\Delta }{\psi }_i,\hfill \end{array}$$(F.38)

$$\begin{array}{cc}\hfill Q{U}_i-U{Q}_i=& P_{i}Psin2\mathrm{\Delta }{\psi }_i.\hfill \end{array}$$(F.39)

We then have

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\delta )=\frac{f_{\mathrm{m}}^2(\delta )}{3P^2}{\displaystyle \sum _{i=1}^N\frac{P_i^2[{sin}^{2}2\mathrm{\Delta }{\psi }_i{tan}^2{\mathrm{\Gamma }}_i+{cos}^{2}2\mathrm{\Delta }{\psi }_i]}{{cos}^2{\mathrm{\Gamma }}_i},}\end{array}$$(F.40)

which can be simplified further, using $P_i^2=p_{\text{max}}^2{I}_i^2{\cos }^4{\Gamma }_i$, to

$$\begin{array}{c}\hfill {\mathcal{S}}^2(\delta )=\frac{f_{\mathrm{m}}^2(\delta )p_{\max }^2}{3P^2}{\displaystyle \sum _{i=1}^N{I}_i^2({sin}^{2}2\mathrm{\Delta }{\psi }_i{sin}^2{\mathrm{\Gamma }}_i+{cos}^{2}2\mathrm{\Delta }{\psi }_i{cos}^2{\mathrm{\Gamma }}_i).}\end{array}$$(F.41)

Furthermore, in our phenomenological model the total intensity is split equally among the N layers, so that I_i = I/N. Therefore,

$$\begin{array}{c}\hfill \mathcal{S}(\delta )=\frac{f_{\mathrm{m}}(\delta )}{\sqrt{3N}}\frac{p_{\max }}{P/I}\mathcal{A},\end{array}$$(F.42)

where 𝒜 is a geometrical factor that depends on the magnetic-field structure in the layers, with

$$\begin{array}{c}\hfill {\mathcal{A}}^2=\frac{1}{N}{\displaystyle \sum _{i=1}^N({sin}^{2}2\mathrm{\Delta }{\psi }_i{sin}^2{\mathrm{\Gamma }}_i+{cos}^{2}2\mathrm{\Delta }{\psi }_i{cos}^2{\mathrm{\Gamma }}_i).}\end{array}$$(F.43)

F.6. Application to Planck data: the case for strong turbulence

In line with our analysis of the data, we compute the mean of 𝒮(δ) over those pixels that have the same polarization fraction p. This gives

$$\begin{array}{c}\hfill {\langle \mathcal{S}(\delta )\rangle }_p=\frac{f_{\mathrm{m}}(\delta )}{\sqrt{3N}}\frac{p_{\max }}{p}{\langle \mathcal{A}\rangle }_p.\end{array}$$(F.44)

Application of the phenomenological model to the Planck data in Planck Collaboration Int. XLIV (2016) shows that good fits are obtained for the parameters p_max = 0.26, f_M = 0.9, and N = 7. This value of f_M implies rather strong turbulence, and therefore the angles Γ_i and Δψ_i are uncorrelated, yielding ${\langle \mathcal{A}\rangle }_{\mathrm{p}}\approx 1/\sqrt{2}$. In that case, the polarization angle dispersion function simply reads

$$\begin{array}{c}\hfill {\langle \mathcal{S}(\delta )\rangle }_{\mathrm{p}}\approx \frac{f_{\mathrm{m}}(\delta )}{\sqrt{6N}}\frac{p_{\max }}{p}.\end{array}$$(F.45)

We reach the important conclusion that the trend 𝒮 ∝ 1/p observed in the Planck data can be reproduced as a generic behaviour that depends only on the statistical properties of the turbulent magnetic field, without invoking changes in properties of the dust or in its alignment with respect to the magnetic field.

We note that the typical value and dispersion of the product 𝒮 × p depend not only on the properties of the turbulence at the scale of the lag, via the $\stackrel{\text{m}}{.}$ parameter, and on the number of layers N, but also on the maximum polarization fraction p_max that the dust can produce. Estimates of the latter are quite sensitive to the uncertainty on the zero level of the total intensity, as discussed in the main text.

For completeness, Fig. F.3 presents the probability density function (PDF) and cumulative density function (CDF) of various distribution functions. The solid curves correspond to the distribution of the mean-normalized product 𝒮 × p/⟨𝒮×p⟩ for our model taken at a resolution of 160′ (models with up to 20′ resolution are similar). Empirically, the corresponding density functions for a Gamma distribution (Hazewinkel 2013) with shape parameter k = 5 and scale parameter θ = 1/5 reproduce these curves well, i.e., this Gamma distribution has similar statistics. These model and empirical density functions are also in reasonable agreement with the PDF and CDF of 𝒮 × p/⟨𝒮×p⟩ for Planck data at the same 160′ resolution.

F.7. Derivation of the expression for f_m(δ)

In our model, each component B_ix, B_iy, and B_iz of the magnetic field vector in layer i is the sum of the uniform field and a realization of a Gaussian random variable on the sphere, with a power-spectrum C_ℓ = C ℓ^α_M, where ℓ is the multipole. As in Planck Collaboration Int. XLIV (2016), we consider that the non-vanishing modes of the turbulent component start at ℓ = 2.

We normalize the turbulent component by imposing that

$$\begin{array}{c}\hfill {\sigma }_B^2=⟨\delta B_{\mathit{ix}}^2⟩+⟨\delta B_{\mathit{iy}}^2⟩+⟨\delta B_{\mathit{iz}}^2⟩=3⟨\delta B_{\mathit{ix}}^2⟩=1,\end{array}$$(F.46)

where δB_ix = B_ix − B_0x, and similarly for the other components. This in turn imposes that for the uniform magnetic field, from the definition of f_M (Eq. (F.4)),

$$\begin{array}{c}\hfill B_0^2=\frac{1}{f_{\mathrm{M}}^2}.\end{array}$$(F.47)

Parseval’s theorem then relates this normalization to the power spectrum by

$$\begin{array}{c}\hfill \frac{1}{2\pi }{\displaystyle \sum _{\ell =2}^{\infty }(\ell +\frac{1}{2})C_{\ell }=⟨\delta B_{\mathit{ix}}^2⟩=\frac{1}{3}.}\end{array}$$(F.48)

The maps of B_ix, B_iy, and B_iz are smoothed to a FWHM resolution $\omega =2\sqrt{2log2}\sigma$, where σ is the standard deviation of the smoothing circular Gaussian beam. This results in smoothed maps denoted B_ix, ω, B_iy, ω, and B_iz, ω. Through the Fourier transform, the total power in the turbulent part of each of these smoothed maps is:

$$\begin{array}{c}\hfill \langle \delta B_{ix,\omega }^2\rangle =\langle \delta B_{iy,\omega }^2\rangle =\langle \delta B_{iz,\omega }^2\rangle =\frac{1}{2\pi }\sum _{\ell =2}^{\infty }(\ell +\frac{1}{2})C_{\ell }{exp}^{-{\sigma }^2{\ell }^2},\end{array}$$(F.49)

where δB_ix, ω = B_ix, ω − B_0x, and similarly for the other components. The loss of power at large ℓ associated with the smoothing is clearly seen in Fig. F.4.

The factor f_m(δ) appearing in the expression of 𝒮 × p (Eq. (F.45)) is by definition (Eq. (F.8)) the typical relative fluctuation of the magnetic field at those scales comprised in the annulus between δ/2 and 3δ/2 (with δ = ω/2), i.e.,

$$\begin{array}{c}\hfill f_{\mathrm{m}}^2(\delta )=\frac{1}{B_0^2}{\displaystyle \sum _{k=x,y,z}{\langle {(\delta B_{ik,\omega }(\mathit{r}+{\mathit{\delta }}^{\mathbf{\prime }})-\delta B_{ik,\omega }(\mathit{r}))}^2\rangle }_{\delta /2\le ||{\mathit{\delta }}^{\mathbf{\prime }}||\le 3\delta /2}.}\end{array}$$(F.50)

Through Parseval’s theorem, and using the Fourier transform of this annulus of mean radius δ = ||δ|| and width δ (Gautier et al. 1992), we find (Fig. F.4)

$$\begin{array}{c}\hfill f_{\mathrm{m}}^2(\delta )=\frac{3}{2\pi }{f}_{\mathrm{M}}^2{\displaystyle \sum _{\ell =2}^{\infty }(\ell +\frac{1}{2})C_{\ell }{exp}^{-{\sigma }^2{\ell }^2}[H^2\left(\ell \frac{\delta }{2}\right)-H^2\left(\ell \frac{3\delta }{2}\right)],}\end{array}$$(F.51)

where H(x) = 2J₁(x)/x, with J₁ the Bessel function of the first kind of order one, and δ expressed in radians. Using Eq. (F.51), we compute $f_{\text{m}}^2(\delta )$ at a resolution ω = 80′ (corresponding to δ = 40′) for α_M = −2.5. We find $f_{\text{m}}^2(\delta )=0.0192f_{\text{m}}{(\delta )}^2$, which corresponds graphically to the hatched region in Fig. F.4.

To determine the dependence of $f_m^2$ on the lag δ, we simplify Eq. (F.51) by considering a unique ℓ ∼ 1/δ and a constant C_ℓ = C_1/δ in the sum, i.e., we replace the smooth, Gaussian, cutoff by a step in ℓ. This simplification, while being numerically incorrect, conserves the scaling of the integral with δ as long as ℓ ≫ 1.

Thus

$$\begin{array}{c}\hfill f_{\mathrm{m}}^2(\delta )\propto f_{\mathrm{M}}^2(\frac{2}{\delta }-\frac{2}{3\delta })(\frac{1}{\delta }+\frac{1}{2})C_{1/\delta }.\\ \hfill \end{array}$$(F.52)

We note that working with resolutions of 160′ and less gives 1/δ >  40. The 1/2 term can therefore be neglected compared to 1/δ, yielding $f_{\text{m}}^2(\delta )\propto f_{\text{m}}{(\delta )}^2{\delta }^{-2-{\alpha }_{\text{M}}}$.

Recalling δ = ω/2 to convert to ω and renormalizing to 160′, this analysis yields the following scaling with resolution:

$$\begin{array}{c}\hfill f_{\mathrm{m}}(\omega )=0.164f_{\mathrm{M}}{\left(\frac{\omega }{{160}^{\prime }}\right)}^{-1-{\alpha }_{\mathrm{M}}/2},\end{array}$$(F.53)

valid as long as α_M does not depart too much from −2.5.

F.8. Beam depolarization

In this section, we estimate the effect of the resolution on the polarization fraction by quantifying the depolarization that occurs within the beam. This is important not only for comparing our results at 80′ and 160′ but also for taking into account the effects of the difference in resolution between the Planck polarization data and the starlight polarization that occurs within a pencil beam.

Following our approach in the previous section, we compute the difference in the total energy of the turbulent component, $f_{\text{p}}^2({\omega }_1,{\omega }_2)$, between two given resolutions ω₁ and ω₂ >  ω₁, for a given line of sight. We have

$$\begin{array}{c}\hfill f_{\mathrm{p}}^2({\omega }_1,{\omega }_2)=\frac{3}{2\pi }{f}_{\mathrm{M}}^2{\displaystyle \sum _{\ell =2}^{\infty }(\ell +\frac{1}{2})C_{\ell }({exp}^{-{\sigma }_1^2{\ell }^2}-{exp}^{-{\sigma }_2^2{\ell }^2}),}\end{array}$$(F.54)

where σ₁ and σ₂ are related to ω₁ and ω₂ by $\omega =2\sqrt{2log2}\sigma$, as already mentioned. We compute $f_{\text{p}}^2({\omega }_1,{\omega }_2)$ for ω₁ = 80′, ω₂ = 2ω₁ = 160′ and α_M = −2.5. We find $f_{\text{p}}^2({80}^{\prime },{160}^{\prime })=0.058f_{\text{m}}{(\delta )}^2$ (see also Fig. F.4). Following the same approach as for Eq. (F.52), this yields the following scaling with the resolution ω:

$$\begin{array}{c}\hfill f_{\mathrm{p}}(\omega ,2\omega )=0.285f_{\mathrm{M}}{\left(\frac{\omega }{{160}^{\prime }}\right)}^{-1-{\alpha }_{\mathrm{M}}/2}.\end{array}$$(F.55)

From Eqs. (F.53) and (F.55), we conclude that the factor $f_{\mathrm{p}}^2(\omega ,2\omega )/{\stackrel{\text{m}}{.}}^2(\omega )$ is independent of f_M and ω, and only depends on α_M. For α_M = −2.5, it is equal to 3.02.

We now study the effect of smoothing on the polarization fraction map. For that we note that Stokes Q and U exhibit a power spectrum similar to that of the turbulent component of the magnetic field at ℓ ≫ 1 (see Appendix F.9). From Parseval’s theorem, we therefore have, for each layer i

$$\begin{array}{c}\hfill {\langle p_i^2\rangle }_{\omega }-{\langle p_i^2\rangle }_{2\omega }=k{p}_{\max }^2{f}_{\mathrm{p}}^2(\omega ,2\omega ),\end{array}$$(F.56)

where the factor $f_{\text{p}}^2(\omega ,2\omega )$ comes from the loss of power in the turbulent component of the field between the two resolutions (Fig. F.4) and k is a constant to be determined numerically. The different random realizations δB_i are independent from one another for ℓ ≥ 2. At small enough spatial scales (ℓ ≫ 1), this ensures that the various realizations of Q_i and U_i are also independent from one another. We therefore have

$$\begin{array}{c}\hfill \langle p_{\omega }^2\rangle -\langle p_{2\omega }^2\rangle =k\frac{p_{\max }^2}{N}{f}_{\mathrm{p}}^2(\omega ,2\omega ).\end{array}$$(F.57)

A comparison with our numerical model shows that k ≃ 1 (Fig. F.5 presents the agreement with the model for p when taking k = 1), and therefore

$$\begin{array}{c}\hfill \langle p_{\omega }^2\rangle -\langle p_{2\omega }^2\rangle =\frac{p_{\max }^2}{N}{f}_{\mathrm{p}}^2(\omega ,2\omega ).\end{array}$$(F.58)

Using Eq. (F.45), this yields

$$\begin{array}{c}\hfill \langle p_{\omega }^2\rangle -\langle p_{2\omega }^2\rangle =6\frac{f_{\mathrm{p}}^2(\omega ,2\omega )}{f_{\mathrm{m}}{(\delta )}^2(\omega )}{\langle \mathcal{S}\times p\rangle }_{\omega }^2.\end{array}$$(F.59)

Equation (F.59) quantifies by how much the polarization fraction decreases when smoothing from resolution ω to 2ω, because on average (only) $\langle p_{2\omega }^2\rangle <\langle p_{\omega }^2\rangle$.

We can now generalize to the case of smoothing data from a finer resolution ω/2ⁿ to a resolution ω. In that case, we can compute the beam depolarization by a chain sum,

$$\begin{array}{c}\hfill \langle p_{\omega /2^n}^2\rangle -\langle p_{\omega }^2\rangle ={\displaystyle \sum _{i=0}^{n-1}[\langle p_{\omega /2^{i+1}}^2\rangle -\langle p_{\omega /2^i}^2\rangle ],}\end{array}$$(F.60)

where all averages are taken over the entire map. This yields

$$\begin{array}{cc}\hfill \langle p_{\omega /2^n}^2\rangle -\langle p_{\omega }^2\rangle & \approx 6\times 3.02{\displaystyle \sum _{i=0}^{n-1}{\langle \mathcal{S}\times p\rangle }_{\omega /2^{i+1}}^2,}\hfill \\ \hfill & \approx 18.1{\langle \mathcal{S}\times p\rangle }_{\omega }^2{\displaystyle \sum _{i=0}^{n-1}{\left(2^{i+1}\right)}^{2+{\alpha }_{\mathrm{M}}},}\hfill \\ \hfill & \approx 18.1\times \frac{1-2^{n(2+{\alpha }_{\mathrm{M}})}}{2^{-2-{\alpha }_{\mathrm{M}}}-1}{\langle \mathcal{S}\times p\rangle }_{\omega }^2.\hfill \end{array}$$(F.61)

In summary, the change in squared polarization fraction from the coarser scale ω to the finer resolution ω/2ⁿ is

$$\begin{array}{c}\hfill \langle p_{\omega /2^n}^2\rangle \approx \langle p_{\omega }^2\rangle +18.1\times \frac{1-2^{n(2+{\alpha }_{\mathrm{M}})}}{2^{-2-{\alpha }_{\mathrm{M}}}-1}{\langle \mathcal{S}\times p\rangle }_{\omega }^2.\end{array}$$(F.62)

In the pencil-beam limit (n = ∞), we have

$$\begin{array}{c}\hfill \delta p_{\mathrm{beam}}^2\equiv \langle p_{\mathrm{pencil}}^2\rangle -\langle p_{\omega }^2\rangle \simeq \frac{18.1}{2^{-2-{\alpha }_{\mathrm{M}}}-1}{\langle \mathcal{S}\times p\rangle }_{\omega }^2\end{array}$$(F.63)

From symmetry arguments, we also have

$$\begin{array}{c}\hfill \delta {(Q/I)}_{\mathrm{beam}}^2\equiv {\langle {[{(Q/I)}_{\mathrm{pencil}}-{(Q/I)}_{\omega }]}^2\rangle }_{\omega }\simeq 0.5\delta p_{\mathrm{beam}}^2\end{array}$$(F.64)

and the same for U/I.

F.9. Comparison of the analytical expressions with numerical results and application to pencil beams

In Fig. F.5, we compare our analytical expressions for the mean 𝒮 × p (Eqs. (F.53) and (F.45)) and the rms of p (Eq. (F.62)) as a function of the resolution, with numerical results directly computed from the smoothed Stokes maps of our simulated model from Planck Collaboration Int. XLIV (2016), i.e., f_M = 0.9, α_M = −2.5, N = 7 and p₀ = 26%. There are two aspects of the comparison for the analytical model, the normalization and the dependence on resolution.

For 𝒮 × p, we observe a slight normalization difference (7%) between the analytical and numerical results. More precisely, at 160′, we find 𝒮 × p = 0${}_{.}^{\circ }$34 for the analytical expression and 0${}_{.}^{\circ }$32 for the simulation. Both are nevertheless very close to the observational value of 0${}_{.}^{\circ }$31. A deviation from the trend predicted by Eq. (F.53) is observed at low resolutions, as expected from our demonstration in Appendix F.7.

Concerning the beam depolarization, the decrease in the polarization fraction with a decreasing resolution (larger ω) is well reproduced over two orders of magnitude in resolution: the example shown corresponds to a mean rms of p of 11% over the full sky. Nevertheless, Vansyngel et al. (2017) already noted a small (approximately 0.1) difference between the index α_M characterizing the power spectrum of the turbulent component of the magnetic field in the simulation, and the index α_EE and α_BB recovered from the analysis of the EE and BB power spectra. This is also what we find here: the scaling of f_m and f_p with ω is actually closer to ω^0.18, which would correspond to α_M = −2.36, when the model is produced with α_M = −2.5. We show this scaling as the solid lines in Fig. F.5 and this is why in the rest of the paper we consider a scaling

$$\begin{array}{c}\hfill f_{\mathrm{m}}(\delta )\propto {\delta }^{0.18}\propto {\omega }^{0.18}.\end{array}$$(F.65)

The case of a pencil beam is of interest. Applying Eq. (F.63) to the highest polarization fraction observed at 160′, p_max(ω = 160′)≈20%, we can estimate the corresponding rms polarization over that scale in pencil-beam data, which corresponds to n = ∞ and for which the pre-factor in Eq. (F.63) is 63.9 for α_M = −2.36. Expressing the observed ⟨𝒮 × p⟩_160′ = 0${}_{.}^{\circ }$31 in radians, we find an rms p_max(ω = 0′) = 20.5%. In that particular case, depolarization is expected to be very small because when p = p_max, the magnetic field is already ordered and within the plane of the sky and therefore little subject to depolarization. However, the effect would be stronger for lines of sight characterized by a low polarization fraction. For example, if p(ω = 160′) = 6.0% at 160′ resolution, Eq. (F.63) gives an rms polarization fraction p(ω = 0′) = 7.4% over that scale for pencil-beam data.

Appendix G: Noise and systematics in the comparison of submillimetre and optical polarization data

For our estimation of the emission-to-extinction polarization ratios R_P/p and R_S/V (Sect. 6), all observations should ideally be done in a pencil beam and probe the full line of sight through the Galactic dust. Unavoidable departures from this ideal situation therefore introduce systematic effects on the quantities appearing in these ratios. In this appendix, we estimate these various systematic effects.

G.1. Beam depolarization

Systematic distortions of the submillimetre polarization signal occur due to the averaging of Stokes Q and U in the telescope beam. This does not happen with the pencil beam of optical observations. In Appendix F.8 we demonstrate that beam depolarization produces a negative bias in p with respect to the pencil-beam value, with a scatter around this biased value (Eq. (F.63)). For α_M = −2.36, we find $\delta p_{\text{beam}}^2\simeq 63.9{\langle \S \times P/I\rangle }_{\omega }^2$ (Appendix F.9). As a consequence, to compare optical and submillimetre polarization data, we compensate for this systematic beam depolarization by multiplying all Planck Stokes parameters and uncertainties (P, p, σ_P, σ_p) by the factor $\sqrt{1+\delta p_{\mathrm{beam}}^2/p^2}$. Because Q and U play a symmetric role in P²  =  Q² + U², the same correction factor is applied to Q, U, Q/I, and U/I.

G.2. Background distortion

As mentioned in Sect. 6.2, the optical polarization degree p_V is potentially biased by the difference in length probed along the line of sight, compared to the polarized emission in the submillimetre. Under the assumption of a uniform ISM, this bias due to the background may be corrected using Eq. (18). Recognizing that the background is not uniform, we estimate that the uncertainty on $p_V^{\infty }$ (the ideal measure from the infinity) is proportional to the amount of reddening behind the star ΔE(B − V)^∞, ⋆ = E(B − V)^∞ − E(B − V)^⋆ and to the rms of the polarization degree per unit reddening in the background, i.e.,

$$\begin{array}{c}\hfill {\sigma }_{p_V}^{\infty }\simeq \mathrm{\Delta }E{(B-V)}^{\infty ,\star }\times \mathrm{rms}{\left(\frac{p_V}{E(B-V)}\right)}_{\mathrm{bkgd}}.\end{array}$$(G.1)

This equation replaces Eq. (19) when the background is not uniform. We discuss the appropriate value of the last factor in Appendix G.4.1.

For the polarization fraction p_V/τ_V, no renormalization is needed as both p_V and τ_V are measured from the star. The systematic uncertainty on p_V/τ_V from background distortion in then simply taken as:

$$\begin{array}{c}\hfill \sigma \left(\frac{p_V}{{\tau }_V}\right)=\frac{{\sigma }_{p_V}^{\infty }}{{\tau }_V}.\end{array}$$(G.2)

These equations also apply respectively to uncertainties on $q_V^{\infty }$ and $u_V^{\infty }$ (Eq. (G.1)), and q_V/τ_V and u_V/τ_V (Eq. (G.2)).

G.3. Uncertainties related to the reddening maps

The uncertainty on E(B − V)^⋆ stems from the uncertainty on the PS1-based reddening data and the uncertainty σ_θ^⋆ on the stellar parallax. We estimate the former by correlating the PS1 total reddening, E(B − V)^∞, with the Planck optical depth at 353 GHz converted to a reddening, E(B − V)_τ. Figure G.1 shows that these quantities are remarkably well correlated, with a Pearson correlation coefficient r = 0.98 and little scatter (around 15% of the running mean on average). Lacking more precise information, we assume equal contributions from instrinsic scatter between E(B − V)^∞ and E(B − V)_τ on the one hand and from noise in E(B − V)^∞ on the other hand, so that we take σ_{E(B − V)}/E(B − V) = 0.1 for the PS1-based estimates. The small departure from linearity in this correlation has been interpreted as evidence for dust evolution in the diffuse ISM (Planck Collaboration XI 2014; Planck Collaboration Int. XXIX 2016), i.e., for the modification of the dust optical properties in its lifecycle through the ISM. This aspect will be investigated in its relation to dust polarization properties in a future paper.

The uncertainty σ_θ^⋆ on the parallax leads to an additional uncertainty on E(B − V)^⋆ that can be estimated roughly by considering the variations of E(B − V)^⋆ when the parallax varies from θ^⋆ − σ_θ^⋆ to θ^⋆ + σ_θ^⋆, i.e.,

$$\begin{array}{c}\hfill {\sigma }_{E(B-V)}^{{\theta }^{\star }}=\frac{E{(B-V)}_{{\theta }^{\star }-{\sigma }_{{\theta }^{\star }}}^{\star }-E{(B-V)}_{{\theta }^{\star }+{\sigma }_{{\theta }^{\star }}}^{\star }}{2},\end{array}$$(G.3)

where $E{(B-V)}^{\star }{}_{{\theta }^{\star }-{\sigma }_{{\theta }^{\star }}}$ and $E{(B-V)}^{\star }{}_{{\theta }^{\star }+{\sigma }_{{\theta }^{\star }}}$ are the reddenings to the star obtained for the altered parallaxes θ^⋆ − σ_θ^⋆ and θ^⋆ + σ_θ^⋆, respectively. Gathering the two sources of uncertainty, the total uncertainty ${\sigma }_{E(B-V)}^{\star }$ on E(B − V)^⋆ is then

$$\begin{array}{c}\hfill {\sigma }_{E(B-V)}^{\star }=\sqrt{{\left[0.1E{(B-V)}^{\star }\right]}^2+{\left[{\sigma }_{E(B-V)}^{{\theta }^{\star }}\right]}^2}.\end{array}$$(G.4)

This uncertainty then propagates to the quantities used in determining R_S/V (see Sect. 6.4), e.g.,

$$\begin{array}{c}\hfill \sigma \left(\frac{p_V}{{\tau }_V}\right)\simeq \frac{1}{{\tau }_V}\sqrt{{\left({\sigma }_{p_V}^{\infty }\right)}^2+{\left(p_V\frac{{\sigma }_{E(B-V)}^{\star }}{E{(B-V)}^{\star }}\right)}^2},\end{array}$$(G.5)

and similarly for q_V/τ_V and u_V/τ_V.

G.4. Polarization angle difference

The unbiased comparison between submillimetre and optical measurements also requires an agreement in polarization angles. We recall that for interstellar polarization of starlight, the direction of the projection of the magnetic field on the plane of the sky (B_POS) can be inferred directly from the polarization angle. For polarization in emission at 353 GHz, a rotation by 90° is required. In Fig. G.2, we compare the direction of B_POS inferred from the two tracers. The length of each line is proportional to the S/N value for the corresponding polarization angle. If the same dust is probed in the optical and the submillimetre, the directions should be identical, and we do find that the agreement is quite remarkable, although not perfect.

To quantify this agreement, we define the difference in orientation angles between the Planck (or submillimetre, “S") and optical (or visual, “V") polarization data as Δψ_S/V  = (ψ+90°)  −  ψ_V. In terms of Stokes parameters this can be written as (Planck Collaboration Int. XXI 2015)

$$\begin{array}{c}\hfill \mathrm{\Delta }{\psi }_{\mathrm{S}/\mathrm{V}}=\frac{1}{2}\mathrm{atan}2[(U{q}_V-Q{u}_V),-(Q{q}_V+U{u}_V)].\end{array}$$(G.6)

Figure G.3 presents the histogram of Δψ_S/V, of which two important aspects should be understood, i.e., its width and any shift relative to the expected centring on zero.

G.4.1. Standard deviation

We attempt to explain the width of the histogram Δψ_S/V (Fig. G.3) by simulating what is expected from noise and systematic effects. We start from Planck data at a resolution of 40′ and starlight data, rotated to assume perfect orthogonality. Then we add fluctuations drawn from Gaussian distributions having the estimated variances for each random and systematic uncertainty considered among the following: noise in the submillimetre and optical, beam depolarization (Appendix G.1), and background distortion (Appendix G.2). Finally, with the resulting Stokes parameters, Q and U at 353 GHz, and q_V and uv in the optical, we compute Δψ_S/V for these simulated data.

The width of the observed histogram is too large to be explained by noise only, as shown by the dashed red curve that results from simulating only effects of the noise in the submillimetre and optical data (Fig. G.3). Part of this discrepancy may be accounted for by the dispersion of polarization angles in the optical within a Planck beam, i.e., at scales that cannot be probed by Planck. To estimate this contribution, we analyse the starlight polarization data alone, ignoring Planck data. We compute the standard deviation σ_V/V of the difference in polarization angles in the optical, for those stars at an angular distance δ smaller than 40′ (one FWHM of the Planck beam). We find σ_V/V = 26°. We note that by its nature this dispersion incorporates twice the variance of optical polarization angles compared to that which enters into the Δψ_S/V histogram. Thus the expected standard deviation from this effect would be about 18°. Three random and systematic effects affecting starlight polarization measures could explain that star-to-star dispersion : noise in the optical, turbulence at scales smaller than the Planck beam (Appendix G.1); and background distortion (Appendix G.2). For these three, we consider Gaussian fluctuations having variances ${\sigma }_{p_V}^2$, $\delta p_{\text{beam}}^2$²⁴ and ${({\sigma }_{p_V}^{\infty })}^2$, respectively, the latter depending on an unknown parameter, rms(p_V/E(B − V))_bkgd. We find that simulations including both optical noise and turbulence within the Planck beam produce a star-to-star dispersion of 13°, well below the observed value of 18°. This points to the need for a contribution from background distortion, and thus a need to estimate rms(p_V/E(B − V))_bkgd. The present data sample provides lower and upper bounds to the rms of p_V/E(B − V), respectively 7% for the rms of p_V/E(B − V) in our sample of 1505 stars, and 13% for the maximum observed polarization fraction (see Sect. 6.6). Assuming a value of rms(p_V/E(B − V))_bkgd = 8% close to the rms over our sample of stars, we obtain a dispersion of 12${}_{.}^{\circ }$6 arising from the background distortion ${\sigma }_{p_V}^{\infty }$. The total dispersion including optical noise, turbulence within the Planck beam, and background distortion is then 16°, still significantly below our estimate of 18°. We note that the individual contributions to the angular dispersion do not quite add up in quadrature to the one obtained when taking into account all sources of errors, because even though the fluctuations of Stokes parameters are Gaussian-distributed, the polarization angles are not (Naghizadeh-Khouei & Clarke 1993). A total dispersion of 18° would require rms(p_V/E(B − V)_bkgd) = 11%, a choice that is probably too extreme given the observational constraints, the uncertainties on our estimated value for σ_V/V, and other physical effects not included in this analysis (see Appendix G.4.2 below). Therefore, to be conservative, we will adopt a value of 8% for the rms of p_V/E(B − V)_bkgd to compute the background distortion, which then by itself contributes a dispersion of 12${}_{.}^{\circ }$6.

Our simulated histogram of Δψ_S/V based on optical and submillimetre noise, beam depolarization at 40′ and background distortion is presented in Fig. G.3 (dashed blue line). It is close in shape and width (standard deviation σ_N + S = 20${}_{.}^{\circ }$8) to the observed histogram and by construction it is centred near zero. The contributions to the standard deviation are 13${}_{.}^{\circ }$9 from submillimetre noise, 8${}_{.}^{\circ }$8 from optical noise, 9${}_{.}^{\circ }$6 from beam depolarization, and 12${}_{.}^{\circ }$6 from background distortion. Here again, we warn against the simple quadratic addition of individual sources of uncertainty. We note that this model does not include a contribution from possible Planck systematics, which can be assessed through the E2E simulations presented in Appendix C. For the lines of sight to stars, we built the histogram of the difference in polarization angles between the input maps at 5′ resolution and those at 40′ that went through the simulation pipeline and therefore include estimates of Planck noise and known systematics. The standard deviation of these histograms is 14° ± 1°, depending on the number of simulations used, a value close to that found from Planck noise alone, underlining the low level of remaining systematics in the Planck data.

This analysis of the standard deviation of Δψ_S/V provides some assurance that the subset of uncertainties needed for the quantification of the emission-to-extinction polarization ratios in Sect. 6.4 have been quantified adequately.

G.4.2. Mean difference

The histogram of Δψ_S/V in the data peaks at −3${}_{.}^{\circ }$1, revealing a systematic angle difference (shift) between polarization angle measurements in emission relative to extinction. Given the large number of lines of sight, this shift cannot be explained by a random process, thus pointing to a systematic effect. Although the shift is small and unimportant for evaluating the polarization ratios discussed in Sect. 6, it is potentially important for other investigations. For example, accurate absolute calibration of the polarization angle is critical for future CMB B-mode experiments to avoid systematic effects that could compromise reaching the precision required (Aumont et al. 2018). We have explored the possible origins of such a shift.

There is some evidence for the possibility of background distortion, a systematic difference arising because in the submillimetre dust is observed along a longer path length than probed in the optical. First, within the full 2388-star sample the shift appears to depend on Galactic longitude. Second, when we apply a more stringent criterion on the reddening ratio E(B − V)^⋆/E(B − V)^∞ (see Sect. 6.3), thus reducing the chance of a significant background contamination, we find that the shift is smaller, −2°, and a longitude dependence is no longer evident.

In the above-mentioned Planck E2E simulations (Appendix G.4.1), the mean and median of the difference in polarization angles are consistently shifted from zero by less than 0${}_{.}^{\circ }$25, which shows that the observed difference is not ascribable to any known systematic in the Planck data. Concerning the calibration of the zero point of the polarization angle in the submillimetre, the uncertainties on the orientation of the HFI PSBs at 100, 143, and 217 GHz are below 0${}_{.}^{\circ }$3 (Planck Collaboration Int. XLVI 2016; Planck Collaboration XI 2020), and it seems likely that the uncertainty at 353 GHz is of the same magnitude. This is consistent with multifrequency measurements of the polarization angle of the bright synchrotron emission of the Crab Nebula from 23 to 353 GHz (Aumont et al. 2018).

The calibration of the zero point in the optical is less clear. For the high-latitude polarization surveys considered here, three highly polarized stars were used (Berdyugin et al. 2014). We have examined the extensive historical record of the polarization angle of these stars and find good evidence that each varies with an rms of typically 1${}_{.}^{\circ }$5, but with excursions as large as 5°. Without confidence in sub-degree accuracy of the optical zero point, the possibility remains of a significant contribution to the shift due to this uncertain calibration. There is not a lot of overlap between the optical polarization compilation of Heiles (2000) and the surveys of Berdyugin et al. (2014) used here. However, by a statistical comparison between the polarization angles of stars in one catalogue with those of stars in the other, binned as a function of angular distance, it is possible to investigate any systematic shift. This analysis can be extended by replacing the optical polarization angle measurements of one or both catalogues with the Planck polarization angle at the catalogue positions. This reveals small systematic shifts up to a few degrees that could arise from different zero point calibrations and/or different path lengths probed. From these investigations, a residual contribution to the shift of order 1° seems possible.

A final possibility is actual decorrelation between the submillimetre and optical polarization, e.g., due to a temperature-weighting effect in the submillimetre coupled with a correlation of variations in the heating of aligned grains (or in the dust properties themselves) along the line of sight with variations of the magnetic field orientation (Tassis & Pavlidou 2015).

We cannot pinpoint a single cause for a shift of the magnitude seen. Instead, it seems that several smaller contributions might have conspired to produce the effect observed.

All Tables

All Figures