Planck 2018 results. XII. Galactic astrophysics using polarized dust emission

We present 353 GHz full-sky maps of the polarization fraction $p$, angle $\psi$, and dispersion of angles $S$ of Galactic dust thermal emission produced from the 2018 release of Planck data. We confirm that the mean and maximum of $p$ decrease with increasing $N_H$. The uncertainty on the maximum polarization fraction, $p_\mathrm{max}=22.0$% at 80 arcmin resolution, is dominated by the uncertainty on the zero level in total intensity. The observed inverse behaviour between $p$ and $S$ is interpreted with models of the polarized sky that include effects from only the topology of the turbulent Galactic magnetic field. Thus, the statistical properties of $p$, $\psi$, and $S$ mostly reflect the structure of the magnetic field. Nevertheless, we search for potential signatures of varying grain alignment and dust properties. First, we analyse the product map $S \times p$, looking for residual trends. While $p$ decreases by a factor of 3--4 between $N_H=10^{20}$ cm$^{-2}$ and $N_H=2\times 10^{22}$ cm$^{-2}$, $S \times p$ decreases by only about 25%, a systematic trend observed in both the diffuse ISM and molecular clouds. Second, we find no systematic trend of $S \times p$ with the dust temperature, even though in the diffuse ISM lines of sight with high $p$ and low $S$ tend to have colder dust. We also compare Planck data with starlight polarization in the visible at high latitudes. The agreement in polarization angles is remarkable. Two polarization emission-to-extinction ratios that characterize dust optical properties depend only weakly on $N_H$ and converge towards the values previously determined for translucent lines of sight. We determine an upper limit for the polarization fraction in extinction of 13%, compatible with the $p_\mathrm{max}$ observed in emission. These results provide strong constraints for models of Galactic dust in diffuse gas.

G Fitting methods for the determination of the emission-to-extinction polarization ratios

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 (RATs; 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 1 full-sky maps, whose sensitivity and sky-coverage largely supersede the previously-available Planck Collaboration: Polarized thermal emission from Galactic dust 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)(2010)(2011)(2012)(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 lowfrequency 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 polarizationsensitive channel of Planck, polarized thermal dust emission is about two orders of magnitude stronger than the polarized cosmic microwave background (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 of a set associated with the 2018 release of data from the Planck mission (Planck Collaboration I 2018), we use all-sky maps of dust polarized emission produced from this third public release of Planck data, hereafter referred to as 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 postprocessing 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 S.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 S× 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 (2018).
The paper is organized as follows.In Sect.2, we present the Planck maps of Stokes parameters that are used in the sub-Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).sequent 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 S to the "polarization gradients" commonly used in polarization studies at lower frequencies.Appendix E gives an interpretation of the inverse relationship between the polarization fraction p and the polarization angle dispersion function S, based on a phenomenological model of magnetized interstellar turbulence.Appendix F 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.Finally, Appendix G discusses the fitting methods used in Sect.6 to estimate ratios of the polarized emission to polarized extinction.

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 (2018) and Planck Collaboration III (2018) 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 (2018), 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 (2018), and so we simply recall a few key properties of this postprocessing 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, which is described in Sect.2.2.
Planck Collaboration: Polarized thermal emission from Galactic dust 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 cosmic microwave background (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. 2 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. 3Therefore, 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.Note that in this case, and to avoid oversampling, the output maps are subsequently downgraded from the original HEALPix4 (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 2018) 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 2018), and subtracting that would add noise unnecessarily.

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 2018), deliberately include a model of the CIB monopole.As a first step towards maps suitable for Galactic science, this needs to be subtracted.GNILC postprocessing 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 −1 , or 452 µK CMB using the unit conversion 287.5 MJy sr −1 K −1
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 Hi 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 Hi, as described in Planck Collaboration VIII (2014) 2018).After the GNILC processing, we apply the same Hi correlation procedure to the output maps of I, in particular finding that a Galactic Hi 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 Hi 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 Hii rather than with neutral hydrogen Hi.The fractional contribution might be most important at low Hi 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 19 cm −2 .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 Hi-associated dust, then this value would need to be added to the Galactic Hi offset.On the other hand, part of any dust emission associated with the WIM is probably correlated with Hi 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 Hi offset, gives a fiducial total Galactic offset of 63 µK CMB (corresponding to 0.0181 MJy sr −1 at 353 GHz), to be added back to the GNILC total intensity map at 353 GHz, after the CIB monopole sub-traction.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 −1 ).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 −1 ) and 103 µK CMB (0.0296 MJy sr −1 ), 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.

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 MJy sr −1 K −1 CMB .The resulting GNILC Stokes maps at 353 GHz and uniform 80 resolution are shown5 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 the Sect.3.2 to estimate the statistical uncertainties on the dust polarization properties.
Note that for studies involving the polarization angle dispersion function S (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.

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

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.

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 where the two-argument function atan(−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 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 P, p, and ψ 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 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 ς 2 < 10 −3 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 an 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).

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. 6A 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 equations B. 2 and B.3 of Planck Collaboration Int. XIX (2015), which are valid at sufficiently high S/N in polarization p/σ p .7 Fig. 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). 8he 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 (2018).These E2E simulations are described in detail in appendix A of Planck Collaboration XI (2018).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.

Polarization angle dispersion function
The polarization angle dispersion function S, 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 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 S, while regions where the polarization patterns are more chaotic exhibit larger values, with S = π/ √ 12 ≈ 52 • when the polarization angles are completely uncorrelated spatially.
A map of S 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 S map over the full sky, based on the GNILC-processed Planck 2018 data release at 353 GHz.Because S 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 , S is still significantly biased (see Sect. 4.1).For this reason, we use maps smoothed to 160 and adopt a correspondingly larger lag δ = 80 . 9This map is shown in the top panel of Fig. 4. We computed the statistical uncertainty σ MC S using the Monte Carlo approach discussed in Sect.3.2, but based on the Stokes maps smoothed to 160 resolution.The map of σ MC S 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). 9When considering the Monte Carlo simulations discussed in the previous subsection, we find that the ratio of the ensemble average map S to the map S computed from the smoothed GNILC data has a mean of 0.90 and a median value of 0.97, with a standard deviation of 0.14.For comparison, when working at 80 resolution and a lag of δ = 40 , these values shift to 0.81, 0.87, and 0.19, respectively, which quantifies the bias that remains when working at 80 resolution.

Relationship of S 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 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 S, 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) (6) which encodes only the angular content of the polarization. 10In Appendix D, we show not only that |∇ψ| is better correlated with S than |∇P| is, but also that this can be demonstrated analytically, with the linear dependence of S 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 S.

Noise and bias in S
An estimate of the variance of S due to noise is (Planck Collaboration Int.XIX 2015; Alina et al. 2016): Just like for p, noise on Stokes parameters Q and U induces a bias on S. 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 π/ √ 12 ≈ 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 This expression works well for S/N on S larger than 3, which we ensure by smoothing the Stokes maps.For notational simplicity, in the rest of this paper, we write S to mean the debiased value S db of the polarization angle dispersion function.

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 S. We then examine the joint distributions of p and total gas column density on the one hand, and of S and p on the other hand.Finally, we look into how one striking feature of these maps, i.e., the inverse relationship between S and p, is well reproduced by relatively simple Gaussian models of the Planck polarized sky.The red solid curve corresponds to the "fiducial" Galactic offset for the total intensity, whereas blue and green correspond to the cases of "low" and "high" offset, respectively.The dashed curves show the mean over the 1000 Monte Carlo histograms, and the envelopes shown as coloured regions span the range of the 1000 histograms.
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 in the total intensity offset.We note, however, that the polarization fractions observed reach at least 20 % for any choice of the total intensity offset, putting strong constraints for dust models.
For comparison, the corresponding coloured dashed 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. 11he 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(p i ), and those for the mean histogram over the 1000 realizations, 12 which we write 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 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 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 +3.5 −1.4 ± 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.
Planck Collaboration: Polarized thermal emission from Galactic dust For completeness, Table 1 also gives the same percentiles for the maps smoothed to 160 resolution, showing a further reduction of the bias h(p i ) − h(p).In that case, we find that the "maximum" dust polarization fraction observed is 21.4 +2.2 −1.2 ± 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 E.8, we quantify the effect of smoothing on p and p max in the framework of the analytical model discussed in Sect.4.3.It is found that smoothing from one resolution to another leads to a decrease of p 2 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 E.9).We conclude that our derivation of p max is affected little by the resolution and much more 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 line of sight (see Appendix E and Fig. E.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 on 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., on 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.

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.• 64 ± 0. • 03.This value is compatible with the earlier measurement in Planck Collaboration Int.XIX (2015).

Polarization angle dispersion function
Finally, the distribution function of S 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 S for randomly oriented polarization.As noted by Alina et al. (2016), if an orientation distribution produces a value of S that is somewhat lower (higher) than this, then the addition of noise tends to make S 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 S for different ranges in polarization fraction (p < 1 %, 1 % < p < 5 %, and p > 5 %) for the GNILC data at of the envelope at 160 resolution is compatible with the largest values found in the map of σ MC S (Fig. 4).Fig. 7 shows that, for the case of an 80 FWHM and lag δ = 40 , at large values of S the DF of the mean 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 S, most of the points have low values, underlining the regularity of the polarization angle on large scales.At 160 resolution and a lag of δ = 80 , the distribution of values in the S map for the data peaks at 1. • 7, with mean and median values of 7. 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 .

Two-dimensional distribution functions
In this section we investigate the 2-dimensional joint distribution function 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.

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.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 20 cm −2 ), the maximum polarization fractions (measured for instance by the 99th percentile curve) for the Monte Carlo simulations are 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 of the maximum p as N H increases.
For completeness, we show in Appendix F 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 of 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 in Fig The sharp downturn of the maximum polarization fraction observed near N H ≈ 10 22 cm −2 corresponds to the strong depolarization occurring on lines of sight that probe high column density structures that are not resolved at 80 .

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 S, 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 S 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.As for Fig. 9, we show the running mean of S in each bin of p for the data. 13We show in the next section that simple analytical models suggest that the relationship is actually S ∝ 1/p.Such a trend is shown in Fig. 10 as the dashed white line.

Relationship to models
All of the properties discussed so far, namely the distribution functions of p, ψ, and S, the decrease of the maximum polarization fraction with increasing column density, and the inverse relationship between S and p, are consistent with the analysis first presented in Planck Collaboration Int.XIX Fig. 11.Same as Fig. 10, but for a phenomenological model of the polarized sky, as described in the text.The dashed white line is the same as in Fig. 10.(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 S 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 S and p was also observed in synthetic polarization maps built from numerical simulations of MHD turbulence (Planck Collaboration Int. XX 2015).
In Appendix E 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 0 , 14 the level of the turbulent magnetic field f M relative to the magnitude of the uniform component, and the spectral index α M of this turbulent component.In Fig. 11 we show the 2-dimensional distribution function of S 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 0 = 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 E.9) S p = 0. • 34/p, which is very close to a fit to the observational trend, S p = 0. • 31/p, overplotted in Figs. 10 and 11.
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 S and p is a generic statistical prop-erty that results primarily from the topology of the magnetic field alone.
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 in 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).

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 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) τ ,15 or total column density of hydrogen N H .The calibrations of such rescalings are uncertain and possibly environment dependent.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, and the relation N H = 5.87 × 10 21 cm −2 × E(B− V) τ (Bohlin et al. 1978;Rachford et al. 2009).It is preferable to use τ converted to N H rather than an estimate of the gas column density derived from Hi 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 and the next, 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 .

Origin of the observed variations of the polarization fraction p
The mutual correlations between p, S, and the column density 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 19 cm −2 < N H < 10 22 cm −2 ) 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 20 cm −2 , the mean value is compatible with a constant, then decreases over the range 0.5-2 × 10 21 cm −2 , and eventually becomes rather flat again.Colouring16 the data with S (on a logarithmic scale) we see from the stratification of the data in Fig. 12 that there is a gradient of S mainly perpendicular to the observed trend of p with respect to N H .This analysis indicates that the decreases of p with S and with N H are mostly independent of each other.
Figure 13 shows how the variations of p and S 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 21 cm −2 ).On the other hand, the opposite is seen in S (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 S, i.e., by the magnetic field structure and the depolarization produced by its variations along the line of sight and within the beam.

Exploring beyond first-order trends using S × p
In Sect.4.3, we concluded that the inverse relationship between S and p is a generic statistical property that results quite sim- ply from the topology of the magnetic field alone, and that a trend S ∝ 1/p, close to that observed, is expected on the basis of simple analytical models (Appendix E).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 S × 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, i.e., the maximum polarization fraction p max .As explained in Appendix E, the product S × 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 S, but also S × p with N H , Galactic latitude b, and Galactic longitude l.Note that throughout this entire analysis lines of sight close to the Galactic plane (|b| < 5 • ) are excluded.As expected, the product S × p has smaller variations with N H , b, and l than exhibited by p and S separately, and the decrease of S × p with N H is systematic, without significant departures.
Away from the Galactic plane, the dependence of S × 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 has disappeared in S × p, confirming our interpretation.However, there are still small variations of S × 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 S × p with Galactic longitude are much less pronounced (of the order of 30 %) than those of p and S independently (a factor 3 or so in each case).

Dedicated study for six molecular regions in the Gould Belt
The radiative torques theory (RATs; 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 2015 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).
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 S 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 S and p in the selected molecular clouds.The various curves show the mean S 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 S 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 S 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 S, is accounted for in the framework of our analytical model (see Appendix E.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 S × p scales as FWHM 0.18 (Eq.E.65).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 full-sky phenomenological model described in Sect.4.3 and smoothed to the same resolution of 40 .Our model, which was able to reproduce the trend S ∝ 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 data, that is due to the change in resolution and is already integrated in our expression for S × p, is also observed in the simulation.
As quantified in Appendix E within the framework of the phenomenological model of Planck Collaboration Int.XLIV (2016), the mean value of the S × p product depends on f M , the ratio between the dispersion of the turbulent component of the field and the mean field strength (see Eq. E.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 Hi 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 radiation field G 0 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 S, i.e., look for a correlation between G 0 and S × p.To this end, we use an estimator G R of the radiation field intensity (Guillet et al. 2018;Fanciullo et al. 2015) that is based on the assumption of thermal equilibrium for dust grains.Under this hypothesis, the dust radiance R, which is the integrated intensity of the dust emission (Planck Collaboration XI 2014), is balanced by the heating of dust by a radiation field with intensity G 0 .This radiation field intensity is estimated using the radiance per unit visual extinc-Fig.14.Two-dimensional histograms with background colours encoding the density of points on a logarithmic scale, showing p (top), S (middle), and S × p (bottom) as a function of the column density N H (left), Galactic latitude b (middle), and Galactic longitude l (right).Black curves show the running means calculated as in Fig. 12, with error bars representing the scatter in each bin.For S, which is on a logarithmic scale, the median trend shown (thin black line) follows the density of points more faithfully than does the mean.tion A V , and in practice, the estimate G R is computed through where E(B − V) R 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 the Pan-STARRS1 colour excess map (Green et al. 2018).
Figure 17 shows the correlation between this estimate G R 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 figure 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 d ≈ 18.5 G 1/5.6 R , corresponding to an average temperature of 18.5 K for a standard radiation field G 0 = 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, which leads to G R overestimating G 0 .This plot therefore demonstrates that the dust temperature is in itself a good proxy for the radiation field intensity in these molecular regions.
According to the RAT theory, we would expect warmer grains to be better aligned, and therefore that S × 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 S × 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 .√ n, where σ is the statistical dispersion in the corresponding bin.The dashed horizontal line in the bottom panel is the mean value of S × p at 160 (Fig. 10), corrected for its dependence on the resolution, as per Eq.E.65.Fig. 16.Mean S as a function of p in selected regions in the Gould Belt for the Planck data (top) and for our phenomenological model (bottom, see text), at a resolution of 40 .The black curve indicates the mean trend averaged over all regions.The dashed line is the fit to the mean S = f (p) trend at 160 (Fig. 10), corrected for its dependence on the resolution, as per Eq.E.65.All bins in p contain the same number of pixels, n ≈ 250.Error bars correspond to the uncertainty in the mean, i.e., σ/ √ n, where σ is the statistical dispersion in the corresponding bin.

Multi-resolution view of the variations of p, S, and S × p across the ISM
The discussion of the variations of p, S, and S × 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 S, 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 S 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 S 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 S 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 the column density threshold above which, at each resolution, the mean S/N for p and S is 3, i.e., above which debiased values of p and S are robust.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 S × 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 21 cm −2 , which is the column density threshold for S 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, S, and S × p.At each resolution ω (expressed in arcmin), using Fig. 20 we determine the minimal column density N min H (ω) above which S is correctly debiased, and define a reference interval I N (ω) = [N min H (ω), 10 × N min H (ω)] (indicated by the horizontal colour bars in Fig. 21).For two consecutive resolutions ω and ω/2, we compute the mean values of p, S, and S× p at both resolutions on the common interval I N (ω) ∩ I N (ω/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 nor-Fig.18. Polarization fraction at 353 GHz (top) and product S × p (bottom) as a function of dust temperature, at a resolution of 40 .The black curve indicates the mean trend averaged over all regions.The dashed horizontal line is the mean value of S × p at 160 (Fig. 10), corrected for its dependence on the lag (Eq.E.65).malized 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 S in concerned. 17he mirrored similarity of each detailed variation in the logarithmic representation of p and S 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 of S by the same factor, demonstrating again that the decrease of p with the column density is mainly a  The residual trend in S × p is small, about a 25 % decrease with column density up to 2 × 10 22 cm −2 .For the case |b| > 10 • , the profile of S×p over this same range of N H is quite flat.For the case including |b| > 2 • , S × p decreases systematically with N H .In our phenomenological model, a systematic decrease of S × p is expected at low Galactic latitudes from an increased number of independent layers N along the line of sight, related to the increased dust opacity and/or length probed (Jones et al. 1992), that is not compensated by the inverse effect of increased S due to an increased distance to the observed dust material (recall that S 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 22 cm −2 .
At slightly higher column densities (N H > 3 × 10 22 cm −2 ), we observe a decrease of p, together with an increase of S and S × p.Such a combination cannot be explained by a decrease in grain alignment, which would not affect S.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, all observed variations of the polarization fraction p with N H are inversely related to those of S, a tracer of the depolarization by the turbulent magnetic field.The multiresolution study of the variation of S × 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 22 cm −2 .

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 of the polarization fraction with increasing column density is widely considered as evidence for a systematic de-crease of 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 all variations observed in the polarization fraction p are mirrored in the dispersion of polarization angles, S, a quantity that is independent of the grain alignment efficiency and of dust optical properties.Quantitatively, the near constancy of S × 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 22 cm −2 .The decrease of 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 22 cm −2 .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 any signature in polarization observables, either in the diffuse ISM or in molecular clouds, which 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 S 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 forthcoming paper dedicated to cold cores.

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.This complements our results already presented in Planck Collaboration Int.XXI (2015), limited at that time to translucent lines of sight (0.15 < E(B − V) τ < 0.8) at low Galactic latitudes (|b| < 30 • ).

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 highlatitude stars (Berdyugin et al. 2001;Berdyugin & Teerikorpi 2001, 2002;Berdyugin et al. 2004Berdyugin et al. , 2014)), we extract data for 2461 lines of sight to stars with measured degree of polarization p V , uncertainty σ p V and position 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 (2075 stars with b > 30 • ), and part of the southern cap (316 stars with b < −59 • ).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: 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 the orthographic projection in Fig. 22, 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 position angle.The agreement between the orientations in the optical and in the submillimetre is quite remarkable.
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 = (ψ 353 + 90 • )−ψ V .In terms of Stokes parameters this can be written as (Planck Collaboration Int.XXI 2015) A histogram of ∆ψ S/V is presented in Fig. 23.The histogram peaks at −3 • and its width corresponds to a standard deviation of 26 • .Neither the width, nor the shift of the histogram with respect to zero, can be explained simply by noise in optical and submillimetre polarization data.To understand the physical origin of the discrepant width, we also plot in Fig. 23 a histogram of the optical-to-optical differences in polarization angles ("V/V," in blue) for those stars at an angular distance δ smaller than 40 (one FWHM of the Planck beam); this thus represents the dispersion of polarization angles in the optical within a Planck beam, i.e., scales that cannot be probed by Planck.This dispersion is dominated by small scale fluctuations of these angles through small-scale turbulence and projection effects; noise does not play a role here because of the high S/N in optical polarization data.We note that by its nature the V/V histogram incorporates twice the variance of optical polarization angles compared to that which enters the S/V histogram.We also plot in Fig. 23 a simulated histogram that corresponds to the differences in orientation angles expected from the noise in Q and U at 353 GHz, assuming a perfect orthogonality between optical and submillimetre polarization, and ignoring the noise in p V that is already included in the V/V histogram.Comparing the variances of the underlying distributions, we find that σ 2 S/V ≈ σ 2 V/V /2 + σ 2 sim , as expected, i.e., the variance of ∆ψ S/V can be explained by the small-scale (< 40 ) dispersion of the polarization angles in the optical within a Planck beam, in combination with the Planck noise.
The slight shift of the peak of the histogram (by −3 • ) reveals a systematic angle difference between these polarization measurements in extinction and emission.Such a difference might be expected, for example, if variations in the heating of aligned grains (or in the dust properties themselves) along the line of sight are correlated with variations of the magnetic field orientation (Tassis & Pavlidou 2015).These possibilities will be studied in a future paper.

Estimates for starlight reddening
To compare 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 et al. 2016;Gaia Collaboration 2018) or from the polarization catalogues when Gaia data are not available.Then we derive an estimate of the reddening towards the star, E(B − V) , by interpolating the Pan-STARRS1 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 of Pan-STARRS1 were separately smoothed to a resolution of 40 and downgraded to N side = 256.We also converted the Pan-STARRS1 reddening to the Johnson E(B − V) scale using the relation Schlafly & Finkbeiner 2011).After removing 70 stars falling outside the region covered by Pan-STARRS1 (mainly in the southern hemisphere), there remain 2388 stars for which we have both reddening and optical polarization data.
The uncertainty on E(B − V) stems from the uncertainty on the Pan-STARRS1 reddening data and the uncertainty σ θ on the stellar parallax.We estimate the former by correlating the Pan-STARRS1 total reddening, E(B − V) ∞ , with the Planck optical depth at 353 GHz converted to a reddening, E(B − V) τ .
Figure 24 shows that these quantities are remarkably well correlated, with a Pearson correlation coefficient r = 0.97 and little scatter (around 15 % of the running mean on average).We therefore assume that the uncertainty of for all lines of sight.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., where E(B − V) θ −σ θ and E(B − V) θ +σ θ are the reddening to the star obtained for the altered parallaxes θ − σ θ and θ + σ θ , respectively.The total uncertainty σ E(B−V) of E(B − V) is then simply (17)

Emission-to-extinction polarization ratios
Following the approach in Planck Collaboration Int.XXI ( 2015) for translucent lines of sight (0.15 < E(B − V) < 0.8), we can now derive the following emission-to-extinction polarization ratios for the diffuse high-latitude Galactic ISM (E(B−V) τ < 0.15, corresponding to column densities N H < 10 21 cm −2 ): Here τ V = A V /1.086 is the optical depth to the star in the optical, is the V-band extinction, and R V = 3.1 is the ratio of total to selective extinction (Fitzpatrick & Massa 2007).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.3.2 and 6.4.

Selection of lines of sight
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 25 presents the histogram of the ratio E(B − V) /E(B − V) ∞ of the reddening towards 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.
If we assumed for simplicity that the ISM along the line of sight is uniform (in density, magnetic-field orientation, and dust properties), then the polarization ratio R P/p = P/p V would artificially increase with decreasing E(B − V) /E(B − V) ∞ .Consequently, by neglecting the presence of background material we would overestimate the polarization ratio R P/p by typically 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: with uncertainty We use similar expressions for q ∞ V and u ∞ V and their uncertainties σ ∞ u V and σ ∞ q V .On the other hand, R S/V , as a ratio of non-dimensional quantities, would be unaffected by this uniform background.However, the ISM is not uniform and as a consequence our estimate of both R S/V and R P/p could be biased by the presence of a background whose properties are different from those in the foreground to the star.We minimize this risk by excluding those lines of sight with an important background, as inferred from the ratio E(B − V) /E(B − V) ∞ shown in Fig. 25.We explicitly choose to keep only stars for which E(B−V) /E(B−V) ∞ > 0.75.
Our final sample contains 1656 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.19], with a mean value of 0.03.

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 , u ∞ V ), and derive R S/V through a joint correlation of (Q/I, U/I) with (q V /τ V , u V /τ V ).In Fig. 26 we present the two correlation scatterplots, that for R P/p on the left, and that for R S/V on the right.
The data in emission and extinction should be treated on an equal footing, neither being in any sense a response to the other.Therefore, estimating the slopes of these correlations warrants using tailored algorithms that are in this sense symmetric, as we discuss in Appendix G.For the fitting, we settle on the bisector version 18 of the Bayesian method of Kelly (2007). 18The symmetric estimator of interest is the slope of the bisector of the y = f (x) and x = f (y) regression lines.
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 when 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), without any bias.Furthermore, we model and quantify in Appendix E.8 the potential bias on P, and therefore on the derived polarization ratios, introduced by the beam depolarization present in Planck data at 40 resolution, while starlight polarization (with its pencil-beam) does not suffer from it.We find that this bias does not exceed 3 %.We conclude that beam depolarization can be safely ignored in our analysis, and at the same time confirm the validity of this approximation in our previous study (Planck Collaboration Int. XXI 2015) carried out at a resolution of 7 .
We find a somewhat larger scatter in the correlation plot (Q/I, U/I) with (q V /τ V , u V /τ V ) in Fig. 26 (right), as quantified by the slightly lower absolute value of the correlation coefficient.The fitted slope R S/V = 4.43 ± 0.03 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.5 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 = 4.2, respectively.

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.This is presented in Fig. 27.The sample is divided into 15 independent bins in N H , the dust optical depth at 353 GHz converted to a column density, each bin containing 110 stars.
For the polarization ratio R P/p at low N H , we observe a roughly 30 % increase with increasing N H , from about 4.0 MJy sr −1 at N H ≈ 6 × 10 19 cm −2 to 5.4 MJy sr −1 at N H > 10 20 cm −2 .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 (see Fig. 27).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. 27.All together, the variations of R P/p and R S/V with the column density show a remarkable trend towards the values determined for translucent lines of sight.
The correlation coefficients are high enough at all column densities to bring confidence in our results.We observe a loss of correlation (lower absolute values of the correlation coefficients) below N H < 10 20 cm −2 , an effect that might result in a underestimation of R P/p and R S/V in this range of column densities.
, yielding an estimate of R P/p .Right: Normalized Stokes parameters (Q/I, U/I) versus (q V /τ V , u V /τ V ), yielding an estimate of R S/V .The slopes of the correlations are obtained using the Bayesian fitting method of Kelly (2007) in its bisector version.Both the Pearson correlation coefficient and the correlation coefficient inferred from the Bayesian method are listed.Fig. 27.Emission-to-extinction polarization ratios.The black curves show the ratios R P/p (left) and R S/V for the fiducial offset in I (right), at 40 resolution, as a function of the column density, N H . Red and purple dashed curves represent the Pearson and Bayesian correlation coefficients, respectively (with the scale on the right axis).The results at a resolution of 80 (squares) and 20 (triangles) are also shown.On the right panel, the upper and lower dashed orange curves represent the trend for the low and high offsets in I, respectively, at the reference resolution of 40 .The blue band shows in each panel the mean value, together with its uncertainty domain, found for the range of column densities considered in Planck Collaboration Int.XXI (2015).

Maximum p V /E(B − V) at low column densities
The maximum starlight polarization degree per unit reddening in the optical, p V /E(B − V) < 9 %, was first estimated by Serkowski et al. (1975).It is often considered as 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 its 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.15.
To study the dependence of the polarization fractions P/I and p V /E(B − V) on column density, we combine data for the 1656 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).As in previous sections, we use the MAS estimator (Plaszczynski et al. 2014) to debias  2015), along with the estimates of maximum polarization.The sample in black is the one in our current study (1656 stars).For each sample, we plot the 99th percentile in p and p V /E(B − V), along with the uncertainty on the derivation of this percentile (see text).The fit from Fosalba et al. (2002), corresponding to p V ∝ E(B − V) −0.2 , is shown for comparison.
the degree of polarization p V in the optical and polarized intensity P in the submillimetre.
Figure 28 shows how these polarization fractions in emission and extinction vary with column density.We overplot the upper limit p V ≤ 9%E(B − V) of Serkowski et al. (1975), the nonlinear fit proposed by Fosalba et al. (2002) for polarization in extinction (Fig. 28, 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 19 ), and a polarization ratio R S/V = 4.5 (Fig. 26), 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. 28 (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 %.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 frac- 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 polariza- 19 The percentile curves in Fig. 28 have uncertainties that are computed in the following way.For each N H bin, the values are sorted and the one closest to the 99th percentile of the distribution is taken as the value for the 99th percentile.The uncertainty interval then spans the range between the value just preceding the 99th percentile value and that just following it.tion 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 20 cm −2 ).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 %, i.e., by about the same factor.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 S× p remains constant over the range of column densities probed here.Following our analysis in Sect.5, we therefore attribute the major part of this decrease of the maximum polarization fraction when the column density increases to the increasing depolarizing effect from the structure of the magnetic field along the line sight, with little room for a systematic decrease of 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.

Conclusions
In this paper, we have analysed the Planck 2018 thermal dust polarization data at 353 GHz, in association with Planck data at other frequencies.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 S, 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 max = 22.0 +3.5 −1.4 ± 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).
We confirmed that the statistical properties of p, ψ, and S 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 S. 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 21 cm −2 ) 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 S, which does not depend on the dust physics.
Based on that analysis, we introduced the product S × 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 S × p exhibits smaller and smoother variations than either p or S 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, S ∝ 1/p.Strikingly, the S 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 S, and therefore by the structure of the magnetic field along the line of sight (Planck Collaboration Int. XX 2015).Considering then the product S × 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 S × p does decrease, but only by about 25 % between N H ≈ 10 20 cm −2 and N H ≈ 2 × 10 22 cm −2 , while the polarization fraction p decreases by a factor 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 (1656 stars selected) allowed us to study the de-pendence of these polarization ratios on column density.We found R P/p to increase rapidly with N H , from 4.2 MJy sr −1 at N H = 6 × 10 19 cm −2 to 5.4 MJy sr −1 for N H > 10 20 cm −2 , the same value as for translucent lines of sight.The polarization ratio R S/V was found to be compatible on average (around 4.5) with that for translucent lines of sight (4.2 ± 0.5), having a decreasing or flat 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.

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 polar-ization 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 decreases 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 S, 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 −1 , 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 into account these constraints.
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 20 cm −2 to 10 22 cm −2 .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, starforming filaments (N H ≈ 10 22 cm −2 ), 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 of 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.XXXV (2016), we studied the relative orientation between filamentary structures of matter and the magnetic field towards molecular clouds of the Gould Belt, Planck Collaboration: Polarized thermal emission from Galactic dust probing lines of sight with total gas column densities N H from around 1021 to 10 23 cm −2 , 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 about N H ≈ 2 × 10 21 cm −2 .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 T E 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 0 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 from 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 0 = 26 ± 3 %), and the ratio between the strengths of the turbulent and mean components of the GMF ( f M = 0.9 ± 0.1).This phenomenological framework was further advanced by Ghosh et al. (2017) and Vansyngel et al. (2017)., 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 : where the sum extends over the seven Planck polarizationsensitive 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. 20The 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 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 halfdifference of the Planck half-mission maps X i,HM 1 and X i,HM 2 : 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 where the 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 N 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 , 21 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 multiresolution 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 also related to the resolution of the Stokes maps by a factor of √ 2. We then compute the local (co)variance value in each region of a given resolution of the Planck and GNILC noise maps, N 353 GHz and N GNILC , for instance: 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 renormalized as 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 K 2 CMB to MJy 2 sr −2 using the conversion factor 287.5 MJy sr −1 K −1 CMB at 353 GHz (Planck Collaboration III 2018).
Figure B.2 shows these covariance maps for the variable resolution case, while Fig. B.3 shows the same maps at the common, uniform 80 resolution.The sky patterns of these uniform resolution covariance maps are extremely similar to those directly taken from the Planck 2018 data, with a significant improvement in the amplitudes, as shown in Table B.1, which gives the means and standard deviations of the ratios σ XY,GNILC /σ XY between the GNILC covariance maps σ XY,GNILC at 80 resolution and N side = 128 to the corresponding maps σ XY in the Planck 2018 data release.The GNILC covariances are overall smaller, and the ratio between two corresponding maps has a very small scatter, underlining the similarity between the patterns in the two maps.Consequently, we use these GNILC-processed covariance maps to assess statistical uncertainties in our analysis.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.ence appearing in the definition of S: 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 on the where the spatial average takes into account that l and θ are independent variables.The former simply yields and the latter average is over θ ∈ [0, 2π].In that average, taking the square and averaging over θ cancels the mixed product, so that 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  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 ) provided that the ratio f m (δ)/ cos Γ i is still small, which only fails if the direction of the mean magnetic field is very close to the line of sight.

E.4. Polarization angle and Stokes parameters
The angle φ i is related to the polarization angle ψ where I i and P i are the total and polarized intensity 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 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 b u 3 C C P o 0 I 5 z q + a 9 e T p x u a z r e f 1 7 Z 3 d v f 3 G i 5 f j I i t l K E Z h F m f y N P A L E U e p G K l I x e I 0 l 8 J P g l h M g u t 3 O j / 5 J G Q R Z e k H t c j F e e J f p t F F F P q K 0 K z x e h p k 8 b x Y J F y W g 9 V s m v j q S i Z L t Z o 1 m k 6 r 5 7 l e 1 7 O d V t v r t 9 0 j B m 6 n 7 7 / w H T + s j 9 Y X 6 6 t 1 8 5 d q 1 S r N K / w 3 r N s / Z 8 i a 4 w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " k N 7 M B 6 s 5 u D b u 3 C C P o 0 I 5 z q + a 9 e T p x u a z r e f 1 7 Z 3 d v f 3 G i 5 f j I i t l K E Z h F m f y N P A L E U e p G K l I x e I 0 l 8 J P g l h M g u t 3 O j / 5 J G Q R Z e k H t c j F e e J f p t F F F P q K 0 K z x e h p k 8 b x Y J F y W g 9 V s m v j q S i Z L t Z o 1 m k 6 r 5 7 l e 1 7 O d V t v r t 9 0 j B m 6 n 7 7 / w H T + s j 9 Y X 6 6 t 1 8 5 d q 1 S r N K / w 3 r N s / Z 8 i a 4 w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " k N 7 M B 6 s 5 u D b u 3 C C P o 0 I 5 z q + a 9 e T p x u a z r e f 1 7 Z 3 d v f 3 G i 5 f j I i t l K E Z h F m f y N P A L E U e p G K l I x e I 0 l 8 J P g l h M g u t 3 O j / 5 J G Q R Z e k H t c j F e e J f p t F F F P q K 0 K z x e h p k 8 b x Y J F y W g 9 V s m v j q S i Z L t Z o 1 m k 6 r 5 7 l e 1 7 O d V t v r t 9 0 j B m 6 n 7 7 / w H T + s j 9 Y X 6 6 t 1 8 5 d q 1 S r N K / w 3 r N s / Z 8 i a 4 w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " k N 7 M B 6 s 5 u D b u 3 C C P o 0 I 5 z q + a 9 e T p x u a z r e f 1 7 Z 3 d v f 3 G i 5 f j I i t l K E Z h F m f y N P A L E U e p G K l I x e I 0 l 8 J P g l h M g u t 3 O j / 5 J G Q R Z e k H t c j F e e J f p t F F F P q K 0 K z x e h p k 8 b x Y J F y W g 9 V s m v j q S i Z L t Z o 1 m k 6 r 5 7 l e 1 7 O d V t v r t 9 0 j B m 6 n 7 7         and the polarized emission is assumed to arise from a small number of layers (here N = 3) in which the total magnetic field B = B 0 + B t is the sum of a uniform field B 0 and an isotropic turbulent field B t that is taken, in each layer, as a different realization of a Gaussian random field in three dimensions.
polarization angle is δψ i = δφ i , and so by inserting Eqs.(E.12) and (E.13) for the fluctuations of the angles we obtain 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.

E.5. Polarization angle dispersion function
The polarization angle dispersion function S 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): . (E.20) Because we are interested in the average behaviour of S, we will ultimately consider the mean of this expression over the position of the central pixel as well.
The distribution function of S (Fig. 7) shows that most pixels have a small dispersion of polarization angles, S 10 • .For is the inclination angle of the magnetic field vector B with respect to the plane of the sky (yz), and φ ∈ [0, 2π] is the angle, counted positively clockwise from the North, between the z axis and the projection of the magnetic field vector onto the plane of the sky.The polarization direction is also in the plane of the sky and perpendicular to that projection, making with the z axis an angle All constructions except B and Γ are in the plane of the sky.these values, it is safe to approximate the arctangent by its argument, so that we may write 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, while for the displaced pixel j we can write (E.27) In the latter expression, the second and third terms are most likely negligible compared to the polarized emission P 2 at the central pixel, especially when averaged over index j, and so they can be ignored.We therefore have 4 S 2 (δ) ≈ U δQ( j) − Q δU( j) 2 j P 4 , (E.28) because P is independent of the pixel j.Appearing in the numerator are the averages δQ( j) 2 j , δU( j) 2 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 g 2  Combining the above expressions, we then have (E.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: (E.39) We then have where A is a geometrical factor that depends on the magneticfield structure in the layers, with sin 2 2∆ψ i sin 2 Γ i + cos 2 2∆ψ i cos 2 Γ i .(E.43) E.6.Application to Planck data: the case for strong turbulence In line with our analysis of the data, we compute the mean of S(δ) over those pixels that have the same polarization fraction p.This gives 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 A p ≈ 1/ We reach the important conclusion that the trend S ∝ 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 S × p depend not only on the properties of the turbulence at the scale of the lag, via the f 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.

E.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 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.E.4), Parseval's theorem then relates this normalization to the power spectrum by The maps of B ix , B iy , and B iz are smoothed to a FWHM resolution ω = 2 2 log 2 σ, 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: The factor f m (δ) appearing in the expression of S × p (Eq. E.45) is by definition (Eq.E.8) the typical relative fluctuation of the magnetic field at those scales comprised in the an-nulus between δ/2 and 3δ/2 (with δ = ω/2), i.e., Regarding the dependence on the lag δ, Eq. (E.51) can be rewritten considering ∼ 1/δ and assuming a constant C = C 1/δ in the sum 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 2 m (δ) ∝ f 2 M δ −2−α M .Recalling δ = ω/2 to convert to ω and renormalizing to 160 , this analysis yields the following scaling with resolution: , (E.53) valid as long as α M does not depart too much from −2.5.

E.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 2 p (ω 1 , ω 2 ), between two given resolutions ω 1 and ω 2 > ω 1 , for a given line of sight.We have . (E.54) We compute f 2 p for ω 1 = 80 , ω 2 = 2ω 1 = 160 and α M = −2.5.We find f 2 p = 0.058 f 2 M (see also Fig. E.3).Following the same approach as for Eq.(E.52), this yields the following scaling with the resolution ω: . (E.55) We now consider a map at a given resolution ω and study the effect on the polarization fraction p of smoothing to twice the resolution, 2ω.Following the approach from Sect.E.4, we define for each line of sight the differences between the values of the Stokes parameters Q and U at resolution ω and those at resolution 2ω: Because we are mainly interested in the effect of the structure of the magnetic field on the smoothed maps, we may assume that the total intensity is uniform at the scale 2ω, in which case we can simplify this expression using p = P/I and P i,2ω = p max I i cos 2 Γ i : .The factor f 2 p (ω, 2ω)/ f 2 m (δ) is independent of ω, but depends on α M .For α M = −2.5, it is equal to 2.14.
We can now generalize to the case of smoothing data from a lower resolution ω/2 n to a resolution ω.In that case, we can compute the beam depolarization by a chain sum: In summary, for data at resolution ω/2 n , the rms p over the coarser scale ω is For S× p, we observe a slight normalization difference (7 %) between the analytical and numerical results.More precisely, at 160 , we find S × p = 0. • 34 for the analytical expression and 0. • 32 for the simulation.Both are nevertheless very close to the observational value of 0. • 31.Concerning the beam depolarization, the decrease of the polarization fraction with a decreasing resolution (larger ω) is reasonably well reproduced: the example shown corresponds to p(ω=160 ) = 10 %, which is the mean p 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 noise covariance matrix to remove the noise bias of the slope.However, like all methods based on debiasing, the BCES method works only when the noise is small and well estimated.Hogg et al. (2010) recommend using Bayesian methods that do not suffer from such limitations.In the Kelly (2007) Bayesian method,28 the parameters of the linear fit are obtained from the posterior distribution with flat priors.We take the mean and standard deviation of the posterior distribution of the slope as our estimate of the slope and its uncertainty.This method is not symmetric, but can be made symmetric by defining its bisector version.It also provides a correlation coefficient that takes into account the measurement errors.Our bisector version of

Fig. 1 .
Fig. 1.From top to bottom: GNILC maps of Stokes I, Q, and U, and polarized intensity P at 353 GHz and uniform 80 resolution in Galactic coordinates, centred on the Galactic centre (GC).The Galactic plane (GP) appears clearly in all maps.Note that the scales for I and P are logarithmic, while those for Q and U are linear.

Fig. 2 .
Fig. 2. Polarization maps for the GNILC-processed data at 353 GHz and uniform 80 resolution: polarization fraction p (top left) and associated statistical uncertainty σ p (top right), polarization angle ψ (bottom left) and associated statistical uncertainty σ ψ (bottom right).The pattern in the σ ψ map arises from the Planck scanning strategy.

PlanckFig. 3 .
Fig. 3. Signal-to-noise ratio (S/N) p/σ p for the polarization fraction in the GNILC-processed data at 353 GHz and uniform 80 resolution.Note that the polar view (bottom) uses a range 1 p/σ p 10 to bring out low S/N regions.

Fig. 4 .
Fig. 4. Top: Polarization angle dispersion function S computed from the GNILC-processed data at 353 GHz and uniform 160 resolution, using a lag δ = 80 .Bottom: Statistical uncertainty σ MC S computed from the Monte Carlo simulations on maps with the same 160 resolution and δ = 80 lag.

4. 1 .Fig. 5 .
Fig. 5. Distribution functions of the polarization fraction p in the GNILC data at 353 GHz and uniform 80 resolution.The red solid curve corresponds to the "fiducial" Galactic offset for the total intensity, whereas blue and green correspond to the cases of "low" and "high" offset, respectively.The dashed curves show the mean over the 1000 Monte Carlo histograms, and the envelopes shown as coloured regions span the range of the 1000 histograms.

Fig. 6 .Fig. 7 .
Fig. 6.Distribution function for the polarization angle ψ in Galactic coordinates for the GNILC data at 353 GHz and uniform 80 resolution.The solid curve shows the histogram of the polarization angles computed directly from the GNILC maps, the dashed curve gives the mean of the 1000 Monte Carlo histograms, and the blue region shows the envelope spanned by the 1000 histograms.

Fig. 9 .
Fig. 9. Two-dimensional histogram showing the joint distribution function of the polarization fraction p from the GNILC data, at 353 GHz and uniform 80 resolution, and the total gas column density N H .This plot uses the fiducial total intensity offset.The black lines show the 5th, 95th, and 99th percentiles of the p distribution in each N H bin, as well as the median p in each N H bin.

Fig. 10 .
Fig. 10.Two-dimensional histogram showing the joint distribution function of S and p at 160 resolution, using a lag δ = 80 .The black curve is the running mean of S as a function of the mean p, in bins of ordered p, with each bin containing the same number of pixels.The error bars represent the standard deviation of S in each bin of p.The dashed white line shows our fit S = 0. • 31/p to this running mean.

Fig. 12 .
Fig.12.Polarization fraction p as a function of the total gas column density, N H , coloured by the polarization angle dispersion function S (on a logarithmic scale).The resolution of the data is 160 , in order to limit the bias in S. The black curve is the running mean of p as a function of the mean N H , in bins of N H that contain the same number of pixels.The top, middle, and bottom running means are calculated for the low, fiducial, and high intensity offsets, respectively.

Fig. 13 .
Fig. 13.Polarization fraction p (top panel), and polarization angle dispersion function S (bottom panel) coloured by the dust temperature, T d , as a function of the column density N H .The resolution of the data is 40 .

Fig. 15 .
Fig. 15.Means from 2-dimensional distributions of polarization properties and column density N H , for selected regions in the Gould Belt, at a resolution of 40 : p (top); S (middle); and S × p (bottom).All bins in N H contain the same number of pixels, approximately 250.Error bars correspond to the uncertainty in the mean, i.e., σ/√ n, where σ is the statistical dispersion in the corresponding bin.The dashed horizontal line in the bottom panel is the mean value of S × p at 160 (Fig.10), corrected for its dependence on the resolution, as per Eq.E.65.

Fig. 17 .
Fig. 17.Correlation between dust temperature T d and our estimate G R for the radiation field intensity, in the selected regions, coloured by N H , and for pixels with N H < 5 × 10 21 cm −2 .The red curve is a prediction for a simple model of dust (see text).

Fig. 19 .
Fig. 19.Mean of p (top), S (middle), and S × p (bottom) as a function of N H , for various resolutions, over the full sky (excluding the Galactic plane, |b| > 5 • ).Dotted lines correspond to trends affected by noise bias.Dashed lines correspond to the uncertainty on the total intensity offset, shown only for 160 data.The background is the density of points at a resolution of 10 .

Fig. 20 .
Fig. 20.Mean S/N of p (top), and S (bottom) as a function of N H , for various resolutions, over the full sky (excluding the Galactic plane, |b| > 5 • ).Error bars correspond to the scatter in each bin.The dashed line indicates the minimal S/N that ensures reliable mean values for debiased quantities.

Fig. 21 .
Fig. 21.Mean S, p, and S × p as a function of N H , combining results from Planck maps at optimal resolutions for all lines of sight above |b| > 5 • .For clarity, S has been raised vertically by a factor of 2. Upper and lower red dashed curves show the corresponding values using the low and high total intensity offsets, respectively.In contrast to other plots, the running means are computed here for bins of equal logarithmic size, which therefore do not contain the same number of pixels.Error bars correspond to the uncertainty in the mean, σ/ √ n, where σ is the statistical dispersion and n is the number of lines of sight in the corresponding bin.Results of the same analysis with different selection criteria on Galactic latitude are shown by dashed (|b| > 10 • ) and dotted (|b| > 2 • ) curves.Horizontal coloured bars indicate for each resolution ω the column density interval I(ω) used in the renormalization procedure (see text).The green band highlights a 25 % decrease of S× p with column density up to 2×10 22 cm −2 .

Fig. 22 .
Fig. 22.Comparison of the orientation of the projection of the magnetic field on the plane of the sky, in orthographic projection with the dust optical depth at 353 GHz as the coloured background, from optical data (top panels) and from Planck data at 353 GHz (bottom panels).The line length is proportional to the S/N on the polarization angle.Northern (left panels) and southern (right panels) Galactic hemispheres are shown, with the Galactic centre situated at the top of each map.

Fig. 23 .
Fig. 23.Histogram of the difference in position angles between Planck-derived angles and optical-polarization-derived angles, ∆ψ S/V = (ψ 353 + 90 • ) − ψ V , with its standard deviation and mean value indicated.The same distribution is plotted for a noise simulation (red dashed curve), and for the difference in position angles between stars with angular distance 0 < δ < FWHM in the optical (blue dashed curve); see text for a description.

Fig. 24 .
Fig. 24.Total reddening observed in the optical, E(B − V) ∞ , as a function of E(B − V) τ , the dust optical depth at 353 GHz converted to a reddening.Each bin of the running mean contains the same number (82) of lines of sight.Error bars represent the standard deviation in each bin.The dashed line corresponds to a one-to-one correlation.

Fig. 25 .
Fig. 25.Histogram of the ratio of the reddening towards the star to the total reddening on the same line of sight, E(B−V) /E(B− V) ∞ , as derived from the Pan-STARRS1 3D cube(Green et al. 2018).The red line indicates the median ratio.In practice we retain lines of sight for which E(B − V) /E(B − V) ∞ > 0.75.

Fig. 26 .
Fig. 26.Correlation between Stokes polarization parameters in emission at 353 GHz and in optical extinction, with the colour in the 2-dimensional histogram representing the density of points.Left: Stokes parameters(Q, U) versus (q ∞ V , u ∞ V ), yielding an estimate of R P/p .Right: Normalized Stokes parameters (Q/I, U/I) versus (q V /τ V , u V /τ V ), yielding an estimate of R S/V .The slopes of the correlations are obtained using the Bayesian fitting method ofKelly (2007) in its bisector version.Both the Pearson correlation coefficient and the correlation coefficient inferred from the Bayesian method are listed.

Fig. 28 .
Fig. 28.Polarization fraction in emission at 353 GHz, p = P/I (left) and normalized polarization fraction in optical extinction, p V /E(B − V) (right), as a function of the column density, N H .The sample in blue contains the 206 translucent lines of sight from Planck Collaboration Int.XXI (2015), along with the estimates of maximum polarization.The sample in black is the one in our current study (1656 stars).For each sample, we plot the 99th percentile in p and p V /E(B − V), along with the uncertainty on the derivation of this percentile (see text).The fit fromFosalba et al. (2002), corresponding to p V ∝ E(B − V) −0.2 , is shown for comparison.

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

Fig. B. 1 .
Fig. B.1.GNILC maps of Stokes I (top left), Q (top right), U (bottom left), and effective FWHM (bottom right) at 353 GHz and varying resolution.The discrete values of the effective FWHM are 5, 7 , 10 , 15 , 20 , 30 , 60 , and 80  .Note that the scale of the I map is logarithmic, while the rest are linear.

2 . 3 .
Fig. C.2. Top: Map of the polarization angle ψ(0) for the input sky of the E2E simulations, at 60 resolution.Middle: Map of the difference between the polarization angle averaged over the 100 realizations of the E2E simulations, ψ , and the input polarization angle ψ (0) , at 60 resolution.Bottom: Distribution function over the sky of the difference between the output and input polarization angle.The solid blue curve is the average of 100 histograms of ψ − ψ (0) from the 100 realizations, while the dashed lines with blue shading between show the ±1 σ dispersion.
Fig. C.4.Histogram of the polarization fractions (top), polarization angles (middle), and polarization angle dispersion functions (bottom).The input data are shown by the black curves and the output of the E2E simulations by the solid blue curves (which are the average of 100 histograms from the 100 simulations), while the dashed lines with blue shading between show the ±1 σ dispersion.

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

2 .
Fig. D.2.Left: Two-dimensional histogram representation of the correlation plot between the angular polarization gradient |∇ψ| and the angular dispersion function S from the GNILC-processed Planck data at 353 GHz and 160 resolution, with a lag of 80 for S. Right: Correlation plot between the polarization gradient |∇P| and the angular dispersion function S. In both plots, the solid black curve shows the mean |∇ψ| or |∇P| in a given bin of S.
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 24 t e x i t s h a 1 _ b a s e 6 4 = " J l C G G 2 Y 5 b 2 3 0 p e P T Q + G p b n y W r Z a I U + 6 x p h c H M e + 4 j m C B H 7 G m 9 G X A O g l n T u g G r O 2 O T l S 8 f c u 4 8O P o Q m Y J 6 4 X O M P I H v u d I o q 4 u 3 T j o i y y k Z V y f X J c r V t X S w 5 w F d g 4 q y E c j L r / g E n 3 E 8 J A i B E M E S T i A A 0 F f F z Y s J M T 1 M C a O E / J 1 n G G C E m l T y m K U 4 R A 7 o n l I u 2 7 O R r R X n k K r P T o l o J + T 0 s Q e a W L K 4 4 T V a a a O p 9 p Z s b 9 5 j 7 W n u l t G q 5 t 7 h c R K 3 B D 7 l 2 6 a + V + d q k V i g G N d g 0 8 1 J Z p R 1 X m 5 S 6 p f R d 3 c / F K V J I e E O I X 7 F O e E P a 2 c v r O p N U L X r t 7 W 0 f E 3 n a l Y t f f y 3 B T v 6 p b U Y P t n O 2 d B 6 6 B q W 1 X 7 7 L B S q + e t L m I H u 9 i n f h 6 h h l M 0 0 C R v j k c 8 4 d k 4 N z L j z r j / T D U K u W Y b 3 4 b x 8 A E m l Z V U < /l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J l C G G 2 Y 5 b 2 3 0 p e P T Q + G p b n y W r Z a I U + 6 x p h c H M e + 4 j m C B H 7 G m 9 G X A O g l n T u g G r O 2 O T l S 8 f c u 4 8O P o Q m Y J 6 4 X O M P I H v u d I o q 4 u 3 T j o i y y k Z V y f X J c r V t X S w 5 w F d g 4 q y E c j L r / g E n 3 E 8 J A i B E M E S T i A A 0 F f F z Y s J M T 1 M C a O E / J 1 n G G C E m l T y m K U 4 R A 7 o n l I u 2 7 O R r R X n k K r P T o l o J + T 0 s Q e a W L K 4 4 T V a a a O p 9 p Z s b 9 5 j 7 W n u l t G q 5 t 7 h c R K 3 B D 7 l 2 6 a + V + d q k V i g G N d g 0 8 1 J Z p R 1 X m 5 S 6 p f R d 3 c / F K V J I e E O I X 7 F O e E P a 2 c v r O p N U L X r t 7 W 0 f E 3 n a l Y t f f y 3 B T v 6 p b U Y P t n O 2 d B 6 6 B q W 1 X 7 7 L B S q + e t L m I H u 9 i n f h 6 h h l M 0 0 C R v j k c 8 4 d k 4 N z L j z r j / T D U K u W Y b 3 4 b x 8 A E m l Z V U < /l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J l C G G 2 Y 5 b 2 3 0 p e P T Q + G p b n y W r Z a I U + 6 x p h c H M e + 4 j m C B H 7 G m 9 G X A O g l n T u g G r O 2 O T l S 8 f c u 4 8O P o Q m Y J 6 4 X O M P I H v u d I o q 4 u 3 T j o i y y k Z V y f X J c r V t X S w 5 w F d g 4 q y E c j L r / g E n 3 E 8 J A i B E M E S T i A A 0 F f F z Y s J M T 1 M C a O E / J 1 n G G C E m l T y m K U 4 R A 7 o n l I u 2 7 O R r R X n k K r P T o l o J + T 0 s Q e a W L K 4 4 T V a a a O p 9 p Z s b 9 5 j 7 W n u l t G q 5 t 7 h c R K 3 B D 7 l 2 6 a + V + d q k V i g G N d g 0 8 1 J Z p R 1 X m 5 S 6 p f R d 3 c / F K V J I e E O I X 7 F O e E P a 2 c v r O p N U L X r t 7 W 0 f E 3 n a l Y t f f y 3 B T v 6 p b U Y P t n O 2 d B 6 6 B q W 1 X 7 7 L B S q + e t L m I H u 9 i n f h 6 h h l M 0 0 C R v j k c 8 4 d k 4 N z L j z r j / T D U K u W Y b 3 4 b x 8 A E m l Z V U < /l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J l C G G 2 Y 5 b 2 3 0 p e P T Q + G p b n y W 0 n b s V P u e y F v C E / 4 v B 0 n 3 A 4 c n 7 e c 0 Y X M t + 5 5 k n p R e C P G M e 8 G 9 j D 0 B p 5 r C 4 J 6 u t 5 x z U P a 3 e Z o S H v p m S q 1 S 6 f 4 9 C e k N H B M m o h 4 C c X y N E P l M + U s 0 X n e E + U p 7 z a m 1 c m 9 A k I F 7 g j 9 S z d j / l c n a x E Y o K p q 8 K i m W C G y O j d 3 y d S r y J s b P 6 o S 5 B A T 0 n b s V P u e y F v C E / 4 v B 0 n 3 A 4 c n 7 e c 0 Y X M t + 5 5 k n p R e C P G M e 8 G 9 j D 0 B p 5 r C 4 J 6 u t 5 x z U P a 3 e Z o S H v p m S q 1 S 6 f 4 9 C e k N H B M m o h 4 C c X y N E P l M + U s 0 X n e E + U p 7 z a m 1 c m 9 A k I F 7 g j 9 S z d j / l c n a x E Y o K p q 8 K i m W C G y O j d 3 y d S r y J s b P 6 o S 5 B A T 0 n b s V P u e y F v C E / 4 v B 0 n 3 A 4 c n 7 e c 0 Y X M t + 5 5 k n p R e C P G M e 8 G 9 j D 0 B p 5 r C 4 J 6 u t 5 x z U P a 3 e Z o S H v p m S q 1 S 6 f 4 9 C e k N H B M m o h 4 C c X y N E P l M + U s 0 X n e E + U p 7 z a m 1 c m 9 A k I F 7 g j 9 S z d j / l c n a x E Y o K p q 8 K i m W C G y O j d 3 y d S r y J s b P 6 o S 5 B A T 0 n b s V P u e y F v C E / 4 v B 0 n 3 A 4 c n 7 e c 0 Y X M t + 5 5 k n p R e C P G M e 8 G 9 j D 0 B p 5 r C 4 J 6 u t 5 x z U P a 3 e Z o S H v p m S q 1 S 6 f 4 9 C e k N H B M m o h 4 C c X y N E P l M + U s 0 X n e E + U p 7 z a m 1 c m 9 A k I F 7 g j 9 S z d j / l c n a x E Y o K p q 8 K i m W C G y O j d 3 y d S r y J s b P 6 o S 5 B A T Fig. E.1.Sketch of the phenomenological model of the dust polarized emission.The observer is represented by the central star, and the polarized emission is assumed to arise from a small number of layers (here N = 3) in which the total magnetic field B = B 0 + B t is the sum of a uniform field B 0 and an isotropic turbulent field B t that is taken, in each layer, as a different realization of a Gaussian random field in three dimensions.

2 .
Fig. E.2.Reference frame for our problem.Γ ∈ [−π/2, π/2] is the inclination angle of the magnetic field vector B with respect to the plane of the sky (yz), and φ ∈ [0, 2π] is the angle, counted positively clockwise from the North, between the z axis and the projection of the magnetic field vector onto the plane of the sky.The polarization direction is also in the plane of the sky and perpendicular to that projection, making with the z axis an angleψ = φ − π/2 [π] ∈ [−π/2, π/2].All constructions except B and Γ are in the plane of the sky.
variables g x , g y , and g z are also uncorrelated from one layer to the next, we have δU(

√ 2 .
In that case, the polarization angle dispersion function simply reads

σ 2 BFig. E. 3 .
Fig. E.3.Power spectrum ( + 1)C /(2π) of a turbulent component of index α M = −2.5, as a function of the multipole .The differential energy lost by smoothing the maps from an initial resolution ω 1 = 80 to ω 2 = 160 is filled in orange, representing a fraction f 2 p (ω 1 = 80 , ω 2 = 160 ) of the original power in the turbulent component (see text in Appendix E.8 and Eq.E.54).Shown as hatched is the turbulent energy implied in the calculation of S at a resolution of ω = 80 with δ = 40 , which is a fraction f 2 m (δ = 40 ) ≡ f 2 m (ω = 80 ) of the original power in the turbulent component (see text in Appendix E.7 and Eq.E.51).Both coloured and hatched regions scale with the resolution ω as ω −2−α M .
where δB ix,ω = B ix,ω − B 0x and similarly for the other components.The loss of power associated with the smoothing is clearly seen in Fig. E.3.

2 i
Fig. E.4.Comparison between numerical results based on smoothed maps of our model of the turbulent magnetic field (diamonds) and the application of our analytical expressions for the decrease of the rms of p (red) by depolarization (Eq.E.64) and the increase of S × p (black) with the resolution (Eqs.E.53 and E.45).The fractional difference is less than 10 % for S × p.
Comparison of the analytical expressions with numerical results and application to pencil-beamsIn Fig.E.4,we compare our analytical expressions for the mean S × p (Eqs.E.53 and E.45) and the rms of p (Eq. E.64) 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 0 = 26 %.Note that there are two aspects of the comparison for the analytical model, the normalization and the dependence on resolution.

1 .
Fig. F.1.Two-dimensional histograms showing the joint distribution function of the polarization fraction p from the GNILC data (at 353 GHz and uniform 80 resolution) and the total gas column density N H .The top plot corresponds to the low total intensity offset, while the bottom plot corresponds to the high total intensity offset.The black lines show the 5th, 95th, and 99th percentiles of the p distribution in each N H bin, as well as the median p in each N H bin.
Planck Collaboration: Polarized thermal emission from Galactic dust 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 of cosmology convention (see Planck Collaboration ES 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

Table 1 .
Statistics from the distribution functions (DFs) for p, as percentages.
a The columns are the following, from left to right: h(p) refers to the DF of the data; h(p i ) refers to the DF of the average p map over the 1000 Monte Carlo realizations; h(p i ) refers to the individual Monte Carlo realizations of the p maps (the values listed give the mean and standard deviations over the 1000 realizations); and h(p i ) refers to the average DF over the 1000 realizations, as shown in Fig.5.
Distribution functions of S at 160 resolution and using a lag of δ = 80 , for different ranges of p, using the fiducial total intensity offset.The distribution function for all points is shown in black and for different ranges of p in separate colours.The distribution functions for the different subsets are scaled to the fractional number of points contained in each range.

Table 2 .
Selected molecular regions in the Gould Belt.