Issue 
A&A
Volume 641, September 2020
Planck 2018 results



Article Number  A12  
Number of page(s)  43  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201833885  
Published online  11 September 2020 
Planck 2018 results
XII. Galactic astrophysics using polarized dust emission
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université ParisDiderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
^{2}
APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{3}
African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg, Cape Town, South Africa
^{4}
Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France
^{5}
Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{6}
Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{7}
CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{8}
CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{9}
Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA
^{10}
California Institute of Technology, Pasadena, CA, USA
^{11}
Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^{12}
Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland
^{13}
Département de Physique, École normale supérieure, PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France
^{14}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{15}
Departamento de Física, Universidad de Oviedo, C/ Federico García Lorca, 18, Oviedo, Spain
^{16}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{17}
Department of Mathematics, University of Stellenbosch, Stellenbosch 7602, South Africa
^{18}
Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{19}
Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa
^{20}
Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{21}
Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, Helsinki, Finland
^{22}
Department of Physics, Princeton University, Princeton NJ, USA
^{23}
Department of Physics, University of California, Santa Barbara, CA, USA
^{24}
Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, IL, USA
^{25}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{26}
Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{27}
Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{28}
Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy
^{29}
Dipartimento di Fisica, Università degli Studi di Trieste, Via A. Valerio 2, Trieste, Italy
^{30}
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{31}
European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{32}
European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{33}
Gran Sasso Science Institute, INFN, Viale F. Crispi 7, 67100 L’Aquila, Italy
^{34}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{35}
INAF – OAS Bologna, Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Area della Ricerca del CNR, Via Gobetti 101, 40129 Bologna, Italy
^{36}
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{37}
INAF – Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy
^{38}
INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
^{39}
INAF/IASF Milano, Via E. Bassini 15, Milano, Italy
^{40}
INFN – CNAF, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{41}
INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{42}
INFN, Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{43}
INFN, Sezione di Milano, Via Celoria 16, Milano, Italy
^{44}
INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{45}
INFN, Sezione di Roma 2, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{46}
Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
^{47}
Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{48}
Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{49}
Institut d’Astrophysique de Paris, CNRS (UMR7095), 98bis boulevard Arago, 75014 Paris, France
^{50}
Institute Lorentz, Leiden University, PO Box 9506, 2300 RA Leiden, The Netherlands
^{51}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{52}
Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{53}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{54}
Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{55}
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{56}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA
^{57}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{58}
Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{59}
Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{60}
Kavli Institute for Particle Astrophysics and Cosmology, Physics and Astrophysics Building, 452 Lomita Mall, Stanford, CA 94305, USA
^{61}
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Chiba 2778583, Japan
^{62}
Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France
^{63}
Laboratoire d’Océanographie Physique et Spatiale (LOPS), Univ. Brest, CNRS, Ifremer, IRD, Brest, France
^{64}
Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{65}
Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{66}
Low Temperature Laboratory, Department of Applied Physics, Aalto University, Espoo 00076 Aalto, Finland
^{67}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{68}
Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{69}
NAOCUKZN Computational Astrophysics Centre (NUCAC), University of KwaZuluNatal, Durban 4000, South Africa
^{70}
National Centre for Nuclear Research, ul. A. Soltana 7, 05400 Otwock, Poland
^{71}
Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00716 Warsaw, Poland
^{72}
Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{73}
SISSA, Astrophysics Sector, Via Bonomea 265, 34136 Trieste, Italy
^{74}
San Diego Supercomputer Center, University of California, 9500 Gilman Drive, La Jolla, San Diego, CA 92093, USA
^{75}
School of Chemistry and Physics, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{76}
School of Physical Sciences, National Institute of Science Education and Research, HBNI, 752050 Jatni, Odissa, India
^{77}
School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK
^{78}
School of Physics and Astronomy, Sun Yatsen University, 2 Daxue Rd, Zhuhai, Tangjia, PR China
^{79}
School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{80}
School of Physics, Indian Institute of Science Education and Research Thiruvananthapuram, Maruthamala PO, Vithura, Thiruvananthapuram, 695551, Kerala, India
^{81}
Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby, BC, Canada
^{82}
Sorbonne Université, Observatoire de Paris, Université PSL, École normale supérieure, CNRS, LERMA, 75005 Paris, France
^{83}
Sorbonne UniversitéUPMC, UMR7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{84}
Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, Moscow 117997, Russia
^{85}
Space Science Data Center – Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy
^{86}
Space Sciences Laboratory, University of California, Berkeley CA, USA
^{87}
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{88}
UPMC Univ. Paris 06, UMR7095, 98bis boulevard Arago, 75014 Paris, France
^{89}
Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex4, France
^{90}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received:
16
July
2018
Accepted:
28
February
2019
Observations of the submillimetre emission from Galactic dust, in both total intensity I and polarization, have received tremendous interest thanks to the Planck fullsky maps. In this paper we make use of such fullsky maps of dust polarized emission produced from the third public release of Planck data. As the basis for expanding on astrophysical studies of the polarized thermal emission from Galactic dust, we present fullsky maps of the dust polarization fraction p, polarization angle ψ, and dispersion function of polarization angles 𝒮. The joint distribution (onepoint statistics) of p and N_{H} confirms that the mean and maximum polarization fractions decrease with increasing N_{H}. The uncertainty on the maximum observed polarization fraction, p_{max} = 22.0_{−1.4}^{+3.5}% at 353 GHz and 80′ resolution, is dominated by the uncertainty on the Galactic emission zero level in total intensity, in particular towards diffuse lines of sight at high Galactic latitudes. Furthermore, the inverse behaviour between p and 𝒮 found earlier is seen to be present at high latitudes. This follows the 𝒮 ∝ p^{−1} relationship expected from models of the polarized sky (including numerical simulations of magnetohydrodynamical turbulence) that include effects from only the topology of the turbulent magnetic field, but otherwise have uniform alignment and dust properties. Thus, the statistical properties of p, ψ, and 𝒮 for the most part reflect the structure of the Galactic magnetic field. Nevertheless, we search for potential signatures of varying grain alignment and dust properties. First, we analyse the product map 𝒮 × p, looking for residual trends. While the polarization fraction p decreases by a factor of 3−4 between N_{H} = 10^{20} cm^{−2} and N_{H} = 2 × 10^{22} cm^{−2}, out of the Galactic plane, this product 𝒮 × p only decreases by about 25%. Because 𝒮 is independent of the grain alignment efficiency, this demonstrates that the systematic decrease in p with N_{H} is determined mostly by the magneticfield structure and not by a drop in grain alignment. This systematic trend is observed both in the diffuse interstellar medium (ISM) and in molecular clouds of the Gould Belt. Second, we look for a dependence of polarization properties on the dust temperature, as we would expect from the radiative alignment torque (RAT) theory. We find no systematic trend of 𝒮 × p with the dust temperature T_{d}, whether in the diffuse ISM or in the molecular clouds of the Gould Belt. In the diffuse ISM, lines of sight with high polarization fraction p and low polarization angle dispersion 𝒮 tend, on the contrary, to have colder dust than lines of sight with low p and high 𝒮. We also compare the Planck thermal dust polarization with starlight polarization data in the visible at high Galactic latitudes. The agreement in polarization angles is remarkable, and is consistent with what we expect from the noise and the observed dispersion of polarization angles in the visible on the scale of the Planck beam. The two polarization emissiontoextinction ratios, R_{P/p} and R_{S/V}, which primarily characterize dust optical properties, have only a weak dependence on the column density, and converge towards the values previously determined for translucent lines of sight. We also determine an upper limit for the polarization fraction in extinction, p_{V}/E(B − V), of 13% at high Galactic latitude, compatible with the polarization fraction p ≈ 20% observed at 353 GHz. Taken together, these results provide strong constraints for models of Galactic dust in diffuse gas.
Key words: polarization / magnetic fields / turbulence / dust / extinction / local insterstellar matter / submillimeter: ISM
© Planck Collaboration 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Interstellar dust grains are heated by absorption of the interstellar radiation field (ISRF), the ambient ultraviolet (UV), visible, and nearinfrared radiation produced by the ensemble of stars in the Galaxy. The grains cool via thermal emission, which is in the farinfrared/submillimetre, as determined by the equilibrium temperature corresponding to a balance between absorbed and emitted power. Thermal emission from the larger grains that dominate the mass in the grain size distribution can be modelled as that of a modified blackbody (MBB) with emissivity ϵ_{ν} = κ_{ν}B_{ν}(T_{d}), where the absorption coefficient κ_{ν} depends on the dust properties (Kruegel 2003). The equilibrium temperature is observed to be of order 20 K (Planck Collaboration XI 2014) for the ISRF found in the bulk of the interstellar medium (ISM).
Starlight polarization, discovered by Hall (1949) and Hiltner (1949), was quickly ascribed to differential extinction by aspherical dust grains with a preferential alignment related to the configuration of the interstellar magnetic field (Davis & Greenstein 1949, 1951). Over the years, a number of theories have been put forward to explain how this alignment occurs and is sustained, despite gas collisions (see the review by Andersson et al. 2015). The mechanism favoured currently involves radiative torques acting on grains subject to anisotropic illumination (RAT; see, e.g., Hoang & Lazarian 2016).
For thermal processes, Kirchhoff’s law states that differential extinction implies differential emission and so the submillimetre thermal emission from dust grains is also polarized, orthogonally to that of extinction. Thus, for dust grains aligned with respect to the Galactic magnetic field (GMF), the observed emission is also partially linearly polarized (Stein 1966; Hildebrand 1988). Because the spin axis of a dust particle is perpendicular to its long axis and alignment is statistically parallel to the local orientation of the magnetic field, the polarization of starlight transmitted through interstellar dust reveals the average orientation of the magnetic field projected on the plane of the sky, whereas the direction of polarized emission is rotated by 90° with respect to the magnetic field.
Observations of this submillimetre emission from Galactic dust, in both total intensity and polarization, have drawn strong attention, thanks to the Planck^{1} fullsky maps, whose sensitivity and skycoverage largely supersede the previouslyavailable data from groundbased, balloonborne (e.g., de Bernardis et al. 1999; Benoît et al. 2004), and space observations (e.g., Gold et al. 2011).
Over the course of four years (2009–2013), Planck surveyed the entire sky in nine frequency bands, from 30 GHz to 857 GHz, providing the best maps to date of the cosmic microwave emission, with unprecedented sensitivity, and angular resolutions varying from 30′ at 30 GHz to 4.′8 at 857 GHz. All but the two highestfrequency 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 highfrequency end. At 353 GHz, which is therefore the highestfrequency polarizationsensitive channel of Planck, polarized thermal dust emission is about two orders of magnitude stronger than the polarized cosmic microwave background (CMB; Planck Collaboration I 2016). It is therefore the channel we use to study this Galactic emission, and several Planck papers have already provided analyses of earlier releases of this data to investigate the link between dust polarization and physical properties of the ISM, most notably the structure of the Galactic magnetic field, properties of dust grains, and interstellar turbulence. In Appendix A, we provide a summary of the main results of these Planck papers, to serve as a useful reference.
In this paper, one in a series associated with the 2018 release of data from the Planck mission (Planck Collaboration I 2020), we use allsky maps of dust polarized emission produced from this third public release of Planck data (hereafter the Planck 2018 data release or PR3) to expand on some of these studies of the polarized thermal emission from Galactic dust. More specifically, our analysis first focuses on a refined statistical analysis of the dust emission’s polarization fraction and polarization angle over the full sky, in the fashion of Planck Collaboration Int. XIX (2015) but based on a 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 𝒮. Here we propose an interpretation of this relationship, showing that it is a generic result of the turbulent nature of interstellar magnetic fields. We therefore further analyse the Planck data by considering the product 𝒮 × p, which allows us to search for deviations from this firstorder 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 componentseparation, i.e., the angular power spectra and spectral energy distributions (SEDs) of the E and B modes, we refer the reader to Planck Collaboration XI (2020).
The paper is organized as follows. In Sect. 2, we present the Planck maps of Stokes parameters that are used in the subsequent analysis. In Sect. 3, we present the fullsky 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 firstorder 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 endtoend (E2E) simulations. Appendix D explains the relationship of the polarization angle dispersion function 𝒮 to the polarization gradients commonly used in polarization studies at lower frequencies. Appendix E provides supplementary figures showing how the behaviour of polarization fraction with total gas column density is affected by the uncertainty on the Galactic zero level. Appendix F provides a demonstration of the inverse relationship between the polarization fraction p and the polarization angle dispersion function 𝒮, based on a phenomenological model of magnetized interstellar turbulence. Finally, Appendix G assesses the noise and systematics that affect the data used in the comparison of visible and submillimetre polarization properties (Sect. 6).
2. Processing Planck maps for Galactic science
The Stokes I, Q, and U maps at 353 GHz that we use in this paper are based on products from the Planck 2018 data release. The processing steps applied to the data to compute the Planck 2018 frequency maps are presented in Planck Collaboration II (2020) and Planck Collaboration III (2020) for the Low Frequency Instrument (LFI) and High Frequency Instrument (HFI), respectively. For HFI, the Q and U products used at 353 GHz make use of the polarizationsensitive bolometers (PSBs) only, ignoring the spiderweb bolometer (SWB) data, as recommended in Planck Collaboration III (2020), while the rest, including I at 353 GHz, make use of the complete data set (PSB+SWB).
For our Galactic science applications, we use maps that result from postprocessing with the Generalized Needlet Internal Linear Combination (GNILC) algorithm, developed by Remazeilles et al. (2011); this filters out the cosmic infrared background (CIB) anisotropies, a key feature for Galactic science. These GNILC maps, derived from the Planck 2018 maps, are presented and characterized in Planck Collaboration IV (2020), and so we simply recall a few key properties of this 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 postprocessing to produce these alternative Stokes maps (ASMs) is also described below.
As another important postprocessing step, we need to establish the desired zero level in the total intensity maps for Galactic dust emission, as described in Sect. 2.2.
2.1. GNILC and ASM postprocessing
GNILC is a waveletbased componentseparation method that makes use of both spectral and spatial information to disentangle multidimensional components of the sky emission. In practice it combines data from the different Planck bands and outputs maps at any desired frequency.
In Planck Collaboration Int. XLVIII (2016), GNILC was applied to Planck 2015 total intensity data, effectively separating Galactic thermal dust emission and CIB anisotropies over the entire sky, while simultaneously filtering out noise and CMB contributions. In regions of low dust column density, it was found that the CIB anisotropies are well above the noise, correlated spatially, and provide a significant contribution to the emission power spectrum. We are interested in polarization properties for Galactic dust emission over the full sky, including highlatitude diffuse lines of sight, for which GNILCprocessing 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 fullwidth at halfmaximum (FWHM) increases from the Galactic plane towards high latitudes. The local resolution depends on the local signaltonuisance ratio, which varies differently for intensity and for E and Bmode polarization^{3}. Therefore, the optimal GNILC resolution should, by design, be different for total intensity and for polarization maps. However, for consistency in the astrophysical study of dust intensity and polarization, where the polarization fraction p = P/I is of interest, we adopt a common resolution by imposing that the variable resolution of the GNILC dust maps should be driven by the more stringent signaltonuisance ratio of the Bmode data. In practice, in the Galactic plane the signaltonoise ratio (S/N) in polarization is sufficiently large to allow for the use of the nominal Planck resolution at 353 GHz, while for high Galactic latitude regions data smoothed to 80′ are required.
The GNILC method is also able to provide Stokes maps at a uniform resolution of 80′ over the entire sky, enabling the analysis of polarization properties over the entire sky at a common resolution. It should be noted that in this case, and to avoid oversampling, the output maps are subsequently downgraded from the original HEALPix^{4} (Górski et al. 2005) resolution N_{side} = 2048 to N_{side} = 128.
The equivalent ASM postprocessing step is to subtract the total intensity CMB SMICA map (Cardoso et al. 2008; Planck Collaboration IV 2020) from the Planck 2018 total intensity map at 353 GHz. No subtraction of CIB anisotropies is performed. Compared to the dust signal at 353 GHz, the CMB polarized signal is small, less than 1% (Planck Collaboration IV 2020), and subtracting that would add noise unnecessarily.
2.2. Zero level for total intensity of Galactic thermal dust emission
We recall that Planck had very little sensitivity to the absolute level of emission and so the zero level of the maps of I must be set using ancillary data. This is of central interest for our study, because for the most diffuse lines of sight it directly impacts polarization fractions through p = P/I.
Planck 2018 HFI frequency maps, as delivered (Planck Collaboration III 2020), deliberately include a model of the CIB monopole. As a first step towards maps suitable for Galactic science, this needs to be subtracted. GNILC 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 given in Planck Collaboration III (2020).
This CIBsubtracted 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 21cm emission of H I at low column densities. Nevertheless, this Galactic offset needs to be refined. A favoured method is again based on a correlation of dust emission with H I, as described in Planck Collaboration VIII (2014), Planck Collaboration XI (2014), Planck Collaboration Int. XLVIII (2016), and Planck Collaboration III (2020). After the GNILC processing, we apply the same H I correlation procedure to the output maps of I, in particular finding that a Galactic H I offset of 36 μK_{CMB} should be added to the 353GHz GNILC total intensity map used for polarization at the uniform 80′ resolution. The statistical error of about 2 μK_{CMB} is small compared to the systematic uncertainties that we now discuss.
Because the dust total intensity versus H I correlation has an upward curvature, the estimates of the offset and slope are dependent on the column density range used for the fit. Furthermore, there is an additional source of uncertainty, related to the possibly significant emission from dust that is in the warm ionized medium (WIM), and therefore associated with H II rather than with neutral hydrogen H I. The fractional contribution might be most important at low H I column densities, i.e., in the diffuse ISM.
To assess the systematic effect related to the WIMassociated 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 H Iassociated dust, then this value would need to be added to the Galactic H I offset. On the other hand, part of any dust emission associated with the WIM is probably correlated with H I as well, and in the extreme case of 100% correlation, there would be no correction due to the WIM.
To account for this effect, we adopt a central value of 27 μK_{CMB} which, when added to the Galactic H I offset, gives a fiducial total Galactic offset of 63 μK_{CMB} (corresponding to 0.0181 MJy sr^{−1} at 353 GHz), to be added back to the GNILC total intensity map at 353 GHz, after the CIB monopole subtraction. This fiducial value will be used in the rest of our analysis. It has an uncertainty that we estimate to be 40 μK_{CMB} (corresponding to ±0.0115 MJy sr^{−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 CIBmonopole subtraction is the same. In this case the fiducial Galactic offset is 68 μK_{CMB}, a value that is, not surprisingly, very close to that for GNILC.
2.3. GNILC Stokes maps
For the 353GHz 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 . The resulting GNILC Stokes maps at 353 GHz and uniform 80′ resolution are shown^{5} 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 GNILCprocessed covariance maps σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU} that are used in Sect. 3.2 to estimate the statistical uncertainties on the dust polarization properties.
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. The scales for I and P are logarithmic, while those for Q and U are linear. 
We note that for studies involving the polarization angle dispersion function 𝒮 (Sect. 3.3), we use Stokes maps and covariance maps that are further smoothed to a 160′ FWHM uniform resolution, and downgraded to N_{side} = 64.
2.4. Alternative Stokes maps (ASMs)
For ASMs, as a final step after converting to astrophysical units, we smooth the Stokes I, Q, and U maps uniformly to 10′, 20′, 40′, 60′, 80′, and 160′, downgrading the HEALPix resolution to N_{side} = 1024, N_{side} = 512, N_{side} = 256, N_{side} = 128, and N_{side} = 64, respectively. The covariance matrix maps σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU} are consistently smoothed from Planck 2018 data to the same resolutions using the procedure described in Appendix A of Planck Collaboration Int. XIX (2015).
3. Fullsky thermal dust polarization maps
In this section, we present the maps of Galactic thermal dust polarization over the full sky, derived from the GNILCprocessed Stokes I, Q, and U maps at uniform 80′ resolution.
3.1. Polarization fraction and angle maps
From the GNILC maps of Stokes parameters I, Q, and U at 353 GHz, we build maps of the polarized intensity P, polarization fraction p, and polarization angle ψ. The convention used for the Stokes parameters in the Planck 2018 data release is to measure polarization angles from the direction of the Galactic north and positively towards Galactic west in accordance with the HEALPix convention used in cosmology (see Planck Collaboration 2018, for further discussion). However, as in Planck Collaboration Int. XIX (2015), we conform here to the IAU convention, polarization angles ψ being counted positively towards Galactic east, and so they are computed simply by changing the sign of Stokes U in the Planck 2018 data. Thus
where the twoargument function atan2(−U, Q) is used in place of atan(−U/Q) to avoid the πambiguity. Conversely, the Stokes parameters can be recovered from the total intensity, the polarization fraction, and the polarization angle via
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 , , 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 noisebias parameter that depends on the geometrical properties of the (assumed Gaussian) 2dimensional distribution of the noise in (Q, U) space, assuming a noisefree total intensity I. From the 353GHz GNILC covariance matrices at the uniform 80′ resolution, we can compute this noisebias 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 a MAS estimate of the polarized intensity, P_{MAS}. For notational simplicity, we hereafter drop the subscript “MAS” and write p to mean p_{MAS} and P to mean P_{MAS}. For the GNILCprocessed 353GHz 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.
Fig. 2. Polarization maps for the GNILCprocessed 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. 
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 largescale patterns at higher latitudes can be broadly interpreted as arising from the projection of the local magnetic field in the Solar neighbourhood (Planck Collaboration Int. XLIV 2016; Alves et al. 2018).
3.2. Estimation of uncertainties
There are several types of uncertainties that need to be taken into account in our analysis of the dust polarization maps.
First, there is statistical noise, whose contribution to the uncertainties can be estimated using the covariance maps of the GNILCprocessed 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^{6}. A set of 1000 simulations was computed; results for a set half this size do not change significantly, confirming that 1000 is sufficient. From these simulations we computed 1000 maps of not only p and ψ, but also other derived quantities, such as the polarization angle dispersion function (Sect. 3.3). As discussed in Sect. 4, these are instrumental in detecting any remaining bias (after using the MAS estimator), by investigating whether statistical properties (e.g., the histogram of p) computed on the GNILC maps shown in Fig. 2, are compatible with the ensemble average of the same properties computed on the Monte Carlo simulations. When, as expected, the quantities in the polarization maps are unbiased, the standard deviations of these maps, and of any derived quantity that we compute using the simulations, yield reliable statistical uncertainties.
Using this approach, we computed polarization fraction and polarization angle uncertainty maps σ_{p} and σ_{ψ} (shown in Fig. 2). These are actually very close to the ones obtained using Eqs. (B.2) and (B.3) of Planck Collaboration Int. XIX (2015), which are valid at sufficiently high S/N in polarization p/σ_{p}^{7}. Figure 3 shows the polarization S/N map for the GNILCprocessed 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)^{8}.
Fig. 3. Signaltonoise ratio (S/N) p/σ_{p} for the polarization fraction in the GNILCprocessed data at 353 GHz and uniform 80′ resolution. The polar view (bottom) uses a range 1 ≤ p/σ_{p} ≤ 10 to bring out low S/N regions. 
The statistical absolute uncertainty on polarization fractions is at most 3%, and the statistical uncertainty on polarization angles is completely negligible, at less than 0.1°. Furthermore, based on the results of Montier et al. (2015a), we are confident that the polarization angle bias is less than 10% of this value. Indeed, at 80′ resolution, 99.9% of the sky pixels have an effective ellipticity below 1.25. This quantity characterizes the asymmetry between the noise distributions on Q and U maps in a rotated reference frame that cancels correlated noise between the two. Montier et al. (2015a) show that in this case the bias on the polarization angle is at most of order 7–8% of the statistical uncertainty σ_{ψ}.
Second, we need to estimate the impact of residual systematics arising from the Planck data processing. This is accomplished via a set of 100 endtoend (E2E) simulations that take a model sky as input, simulate the timelines of the instrument taking into account all known systematics, and process these simulated timelines with the mapmaking pipeline described in Planck Collaboration III (2020). These E2E simulations are described in detail in Appendix A of Planck Collaboration XI (2020). The statistical comparison between the input and output polarization maps, which we discuss in Appendix C, shows that the absolute uncertainties from residual systematics are estimated to be ±0.5% on p and ±8° on ψ.
We note that these E2E simulations include realizations of random data noise and so already include part of the statistical uncertainty that is addressed by the Monte Carlo simulations based on the covariance matrices.
Finally, as already mentioned in Sect. 2.2, the quantitative analysis of p towards diffuse lines of sight depends strongly on the value of the Galactic offset used to set the zero level of total intensity for Galactic dust emission. To take this source of uncertainty into account, following the discussion in Sect. 2.2 we consider a fiducial case in which the Galactic offset is 63 μK_{CMB} and also consider a range of ±40 μK_{CMB} about this central value.
3.3. Polarization angle dispersion function
The polarization angle dispersion function 𝒮, introduced in Planck Collaboration Int. XIX (2015) quantifies the local (non)uniformity of the polarization angle patterns on the sky by means of the local variance of the polarization angle map at a certain scale parameterized by a lag δ. It is defined as
where the sum extends over the N pixels, indexed by i and located at positions r + δ_{i}, within an annulus centred on r and having inner and outer radii δ/2 and 3δ/2, respectively. Regions where the polarization angle tends to be uniform exhibit low values of 𝒮, while regions where the polarization patterns are more chaotic exhibit larger values, with when the polarization angles are completely uncorrelated spatially.
A map of 𝒮 at 60′ resolution and using a lag of 30′, based on Planck 2013 data, was shown over a restricted region of the sky in Planck Collaboration Int. XIX (2015). We can now present the 𝒮 map over the full sky, based on the GNILCprocessed Planck 2018 data release at 353 GHz. Because 𝒮 is built from the polarization angle ψ, it is independent of the value chosen for the total intensity offset. However, when computed at uniform 80′ resolution and using a lag δ = 40′, 𝒮 is still significantly biased (see Sect. 4.1). For this reason, we use maps smoothed to 160′ and adopt a correspondingly larger lag δ = 80′^{9}. This map is shown in the top panel of Fig. 4. We computed the statistical uncertainty using the Monte Carlo approach discussed in Sect. 3.2, but based on the Stokes maps smoothed to 160′ resolution. The map of 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).
Fig. 4. Top: polarization angle dispersion function 𝒮 computed from the GNILCprocessed data at 353 GHz and uniform 160′ resolution, using a lag δ = 80′. Bottom: statistical uncertainty computed from the Monte Carlo simulations on maps with the same 160′ resolution and δ = 80′ lag. 
3.4. Relationship of 𝒮 to alternative estimators
Synchrotron studies in the radio domain frequently use another estimator of the uniformity of polarization patterns, the polarization gradient introduced by Gaensler et al. (2011) and defined as
where y and z refer to an orthogonal coordinate system on the plane of the sky. We show in Appendix D that, as far as the Planck thermal dust polarization data are concerned, ∇P is strongly correlated with 𝒮, though not perfectly because of the contribution from the polarized intensity in ∇P. This can be mitigated by considering an angular version of the polarization gradient defined as (Burkhart et al. 2012)
which encodes only the angular content of the polarization^{10}. In Appendix D, we show not only that ∇ψ is better correlated with 𝒮 than ∇P is, but also that this can be demonstrated analytically, with
the linear dependence of 𝒮 on the lag being revealed simply through a firstorder Taylor expansion. We do not use this estimator ∇ψ in the rest of this paper, but note that in practice it might be easier to compute than 𝒮.
3.5. Noise and bias in 𝒮
An estimate of the variance of 𝒮 due to noise is (Planck Collaboration Int. XIX 2015; Alina et al. 2016):
Just like for p, noise on Stokes parameters Q and U induces a bias on 𝒮. Unlike for p, however, this bias can be positive or negative, depending on whether the true value is, respectively, smaller or larger than the value 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 𝒮 larger than 3, which we ensure by smoothing the Stokes maps. For notational simplicity, in the rest of this paper, we write 𝒮 to mean the debiased value 𝒮_{db} of the polarization angle dispersion function.
4. Statistics of thermal dust polarization maps
In this section, we provide a statistical analysis of the quantities represented in the Galactic thermal dust polarization maps derived above. We start by discussing the distribution functions of p, ψ, and 𝒮. We then examine the joint distributions of p and total gas column density on the one hand, and of 𝒮 and p on the other hand. Finally, we look into how one striking feature of these maps, i.e., the inverse relationship between 𝒮 and p, is well reproduced by relatively simple Gaussian models of the Planck polarized sky.
4.1. Distribution functions for p, ψ, and 𝒮
4.1.1. Polarization fraction
The distribution function (DF, or histogram) for p over the full sky is shown in Fig. 5 The solid red curve is the histogram for the GNILC map of p for the fiducial offset in I, while the solid blue and green curves are the corresponding histograms for the low and high offsets, respectively. These clearly show the significant effect induced by the uncertainty on the total intensity offset. We note, however, that the polarization fractions observed reach at least 20% for any choice of the total intensity offset, setting strong constraints for dust models.
Fig. 5. Distribution functions of the polarization fraction p in the GNILC data at 353 GHz and uniform 80′ resolution. The solid red 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. 
For comparison, the corresponding dashed coloured curves are the means of the DFs from the 1000 Monte Carlo simulations. Compared to the solid curves, there is only a small bias, shifting the DF towards higher p in the tail of the distribution. This is less pronounced for the green curves (high total intensity offset) because for this case the statistical changes in I are less important^{11}.
The coloured regions encompassing the mean histograms show the minimum and maximum values of the histogram for any given bin of p over the 1000 samples, i.e., the envelope within which all 1000 histograms lie. Lines defining the edges of the envelope would themselves not be distribution functions; however, they give an idea of the possible spread of the p histograms with varying noise realizations.
It is of interest, for dust models in particular, to estimate p_{max}, the largest value of p over the full sky. To estimate this and provide further quantification, we compute, for each of the total intensity offset values, the 90th, 95th, 98th, 99th, and 99.9th percentiles for the GNILC histogram from the data, which we write as h(p), and for each of the 1000 Monte Carlo histograms, which we write as h(p_{i}), with 1 ≤ i ≤ 1000. From the latter we calculate the mean and the standard deviation, which gives an estimate of the statistical uncertainty of p_{max} in a single realization, such as the data.
These numbers are given in Table 1, alongside the corresponding values for the average p map over the 1000 Monte Carlo realizations, which we write , and those for the mean histogram over the 1000 realizations^{12}, which we write . 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 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 . 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%.
Statistics from the distribution functions of p, given as percentages.
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 for the maximum dust polarization fraction observed at 80′ resolution and 353 GHz over the full sky, where the first uncertainty relates to the systematic effect of the total intensity offset and the effects of residual systematics, and the second covers the statistical uncertainty estimated from the 1000 Monte Carlo realizations.
For completeness, Table 1 also gives the same percentiles for the maps smoothed to 160′ resolution, showing a further reduction of the bias . In that case, we find that the maximum dust polarization fraction observed is . This value of p_{max} and the debiased value at 80′ agree quite well. This shows that smoothing has little effect on the polarization fraction. Of course, the amount of smoothing applied should not be excessive, because of the potential impact of beam depolarization at higher FWHM. In Appendix F.8, we quantify the effect of smoothing on p and p_{max} in the framework of the analytical model presented in Sect. 4.3. It is found that smoothing from one resolution to another leads to a decrease in p^{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 F.9). We conclude that our derivation of p_{max} is not so much affected by the resolution and much more so by the offset in I.
These results are consistent with the finding of Planck Collaboration Int. XIX (2015) that p_{max} > 19.8% at 60′ resolution over a smaller fraction of the sky. We have also checked that they are not significantly affected when selecting only those pixels on the sky for which the S/N in polarization is p/σ_{p} > 3.
As was pointed out in Planck Collaboration Int. XIX (2015) and Planck Collaboration Int. XX (2015), the level of observed polarization fractions is strongly dependent on the angle Γ of the mean magnetic field with respect to the plane of the sky (see Appendix F and Fig. F.2). The distribution function of p must depend on this mean orientation of the Galactic magnetic field. Compared to what would be obtained for a mean field that is everywhere in the plane of the sky, the distribution should be more peaked towards lower values, as we do observe, but the value of p_{max} might still be high, reflecting those parts of the sky with a favourable orientation, i.e., in the plane of the sky. Although the estimate of p_{max} based on percentiles would be impacted, such an analysis (requiring a model of the largescale GMF) is beyond the scope of this paper.
4.1.2. Polarization angle
Figure 6 shows the distribution function for the polarization angle ψ, for which the value of the total intensity offset is unimportant. The comparison between the histogram for the GNILC map of ψ and the mean histogram over the Monte Carlo realizations shows that there is virtually no noise bias. The histograms peak around 0°, which corresponds to an orientation of the GMF parallel to the Galactic plane. Quantitatively, over the 1000 Monte Carlo samples, the ensemble average of the mean polarization angle is −064 ± 003. This value is compatible with the earlier measurement in Planck Collaboration Int. XIX (2015, see their Fig. 3).
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 data, the dashed curve gives the mean of the 1000 Monte Carlo histograms, and the blue region shows the envelope spanned by the 1000 histograms. 
4.1.3. Polarization angle dispersion function
Finally, the distribution function of 𝒮 is shown in Fig. 7. Results for the case of a 160′ FWHM and lag δ = 80′ are shown in green, and for the case of a 80′ FWHM and lag δ = 40′ in blue. As above for p and for ψ, the solid lines are for the GNILC maps, the dashed lines are the Monte Carlo means, and the coloured regions show the span of histograms for the 1000 Monte Carlo realizations.
Fig. 7. Distribution functions of the polarization angle dispersion function 𝒮 in the GNILC data at 353 GHz. The cases shown are for the 160′ resolution using a lag δ = 80′ (in green), and for the 80′ resolution using a lag δ = 40′ (in blue). The solid curves show the histograms computed directly from the GNILC maps, the dashed curves give the mean histogram from the 1000 Monte Carlo realizations for each case, and the coloured regions show the envelope. The dashed vertical line indicates the value corresponding to a completely random polarization pattern. 
It is interesting that these distributions have a tail passing through 52°, the value of 𝒮 for randomly oriented polarization. As noted by Alina et al. (2016), if an orientation distribution produces a value of 𝒮 that is somewhat lower (higher) than this, then the addition of noise tends to make 𝒮 larger (smaller), towards 52°. This tail in the full DF in Fig. 7 is strongly associated with regions where p is small and more susceptible to the influence of noise, as is apparent in Fig. 8, which shows the distribution function of 𝒮 for different ranges in polarization fraction (p < 1%, 1%< p < 5%, and p > 5%) for the GNILC data at 160′ resolution and with a lag δ = 80′. The large values of 𝒮 are also associated with large values of the scatter , as shown by the widening of the envelope at high values of 𝒮 in Fig. 7. The width of the envelope at 160′ resolution is compatible with the largest values found in the map of (Fig. 4).
Fig. 8. Distribution functions of 𝒮 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. 
Figure 7 shows that, for the case of an 80′ FWHM and lag δ = 40′, at large values of 𝒮 the mean DF of the Monte Carlo realizations is clearly biased with respect to the distribution function of the data, which does not even fit within the region spanned by the 1000 Monte Carlo simulations. On the other hand, for 160′ FWHM and lag δ = 80′, the bias is much less apparent and so we focus on the results for this case. Despite the long tail at large 𝒮, most of the points in this tail have low occurrence rates, underlining the regularity of the polarization angle on large scales. At 160′ resolution and a lag of δ = 80′, the distribution of values in the 𝒮 map for the data peaks at 17, with mean and median values of 76 and 46, respectively. The same characteristic values over the 1000 Monte Carlo simulations are, respectively, 19 ± 06, 829 ± 001, and 512 ± 001. Using the 99th percentile, most of the points in the data have 𝒮 ≤ 436, while the Monte Carlo simulations give an estimate of 453 ± 02. We give these values for reference in the future, for instance in work comparing Planck data with MHD simulations and analytical models.
We stress that while the smoothing to 160′ is warranted here for studies including the highlatitude sky, this requirement for smoothing should not be generalized. Indeed, when the analysis is restricted to the approximately 42% of the sky considered in Planck Collaboration Int. XIX (2015), we find that no such bias exists when working at 80′ FWHM and lag δ = 40′. Incidentally, this confirms the results shown in Planck Collaboration Int. XIX (2015) at 60′ resolution and δ = 30′.
4.2. Twodimensional distribution functions
In this section we investigate the 2dimensional joint distribution functions of polarization fraction p and another variable. Therefore, instead of simply presenting a scatter plot, we display a 2dimensional histogram made by binning in the two dimensions and encoding the number in each bin by colour.
4.2.1. Polarization fraction versus total gas column density
In Fig. 9 we display the 2dimensional histogram of p and total gas column density N_{H}, using the GNILC data, at 353 GHz and uniform 80′ resolution, with the fiducial total intensity offset, over the full sky. We discuss the determination of N_{H} at the beginning of Sect. 5. No cut in either S/N or Galactic latitude has been performed here. The colour scale encodes the logarithm of counts in each bin, while the black curves show the 5th, 95th, and 99th percentiles of the p distribution in each N_{H} bin, as well as the median polarization fraction in each N_{H} bin.
Fig. 9. Twodimensional 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. 
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 2dimensional 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 slightly higher than the corresponding values from the data. As expected, this bias is in the same sense as discussed in Sect. 4.1.1 for the distribution function of p. Recall that this is for 80′ resolution; when working at 160′ resolution this bias disappears.
The joint (N_{H},p) distribution has qualitatively the same behaviour as that found over a smaller fraction of the sky in Planck Collaboration Int. XIX (2015): a large scatter of p towards diffuse lines of sight and a decrease in the maximum p as N_{H} increases.
For completeness, we show in Appendix E the effect of the total intensity offset. It is negligible at the high intensity end, where the histograms are similar whether we use the fiducial, high, or low offset values. At the low intensity end, on the other hand, the effect of the offset is more marked. There is a significant increase in characteristic values (highest percentiles) of p for decreasing N_{H} when taking the low offset, and conversely a marked decrease in the maximum p when taking the high offset.
One might wonder if it would be possible to constrain the offset by assuming that p should reflect dust properties at low column densities, and therefore that the offset should be such that the maximum p is approximately constant at low N_{H}. In this respect, the fiducial offset seems more adequate than either the high or low cases, as can be seen by comparing Fig. 9 with Fig. E.1.
The sharp downturn of the maximum polarization fraction observed near N_{H} ≈ 10^{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′.
4.2.2. Polarization angle dispersion versus polarization fraction
In Planck Collaboration Int. XIX (2015), we discovered an inverse relationship between the polarization fraction p and the polarization angle dispersion function 𝒮, working with data over approximately 42% of the sky, at a resolution of 60′ and a lag of δ = 30′. We have verified quantitatively on the same sky region and using the same methodology that the same inverse relationship holds with the Planck 2018 data release; the maps of polarization are very similar where the S/N is high, as expected. In this limited sky region, we also find that the results are only slightly dependent on the adopted Galactic offset.
Extending to the full sky at 160′ resolution and a lag of δ = 80′, we present the 2dimensional histogram of the joint distribution function of 𝒮 and p in Fig. 10. The data clearly show that the inverse relationship seen at low and intermediate Galactic latitudes in Planck Collaboration Int. XIX (2015) persists in the highlatitude sky. In Fig. 10 we also display the running mean of 𝒮 in each bin of p for the data^{13}. We show in the next section that simple analytical models suggest that the relationship is indeed ⟨𝒮⟩_{p} ∝ 1/p. Such a trend is shown in Fig. 10 as the dashed white line.
Fig. 10. Twodimensional histogram showing the joint distribution function of 𝒮 and p at 160′ resolution, using a lag δ = 80′. The black curve is the running mean of 𝒮 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 𝒮 in each bin of p. The dashed white line shows our fit 𝒮 = 031/p to this running mean. 
4.3. Relationship to models
All of the properties discussed so far, namely the distribution functions of p, ψ, and 𝒮, the decrease in the maximum polarization fraction with increasing column density, and the inverse relationship between 𝒮 and p, are consistent with the analysis first presented in Planck Collaboration Int. XIX (2015). Subsequently, phenomenological models of the polarized sky incorporating Gaussian fluctuations of the magnetic field have been developed (Planck Collaboration Int. XLIV 2016; Ghosh et al. 2017; Vansyngel et al. 2017; Levrier et al. 2018). Interestingly, although these models were built to reproduce some 1 and 2point statistics of polarized emission maps, they were not tailored to reproduce the inverse relationship between 𝒮 and p that is evident in the Planck data, and yet they are able to do so very readily and robustly. A similar inverse relationship between 𝒮 and p was also observed in synthetic polarization maps built from numerical simulations of MHD turbulence (Planck Collaboration Int. XX 2015).
In Appendix F we present an analytical derivation of this property, based on the most basic version of these phenomenological models. In that framework, the emission is assumed to arise from a small number N of layers, each emitting a fraction 1/N of the total intensity, and harbouring a magnetic field that is the sum of a uniform component and a turbulent Gaussian component. The main parameters of the model, besides N, are the intrinsic polarization fraction p_{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 the spatial distribution of this turbulent component.
In Fig. 11 we show the 2dimensional distribution function of 𝒮 and p at 160′ resolution, using a lag δ = 80′, for a polarized sky built from such a Gaussian phenomenological model, with α_{M} = −2.5, f_{M} = 0.9, N = 7, and p_{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 F.9) ⟨𝒮⟩_{p} = 034/p, which is very close to a fit to the observational trend, ⟨𝒮⟩_{p} = 031/p, overplotted in Figs. 10 and 11. We show in Appendix F.9 that this relationship depends weakly on the resolution ω according to
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. 
Because changes of dust properties or dust alignment are not included in these phenomenological models nor in the synthetic observations from MHD simulations, we conclude that the inverse relationship between 𝒮 and p is a generic statistical property that results primarily from the topology of the magnetic field.
We also note that neither the phenomenological model of Planck Collaboration Int. XLIV (2016), nor the MHD simulation in Planck Collaboration Int. XX (2015), account for the 3D structure of the ordered (mean) component of the GMF on large scales. The imprint of this ordered component on the dust polarization can be readily identified in the map of the dust polarization angle (Fig. 2). It also impacts the polarization fraction map on large angular scales and thereby the dependency of p on Galactic latitude. For synchrotron polarization, this has been quantified by Page et al. (2007) and MivilleDeschênes et al. (2008) using Galaxywide models of the GMF. As discussed in Alves et al. (2018), a comprehensive model of dust polarization would also need to take into account the structure of the GMF on the scale of the Local Bubble (100–200 pc).
5. Insight from interrelationships and Galactic context
Further insight into statistical measures of the polarization can be gained not only by considering them in relation to one another, but also by studying how they jointly vary with other physical parameters such as dust temperature T_{d} or column density, and how these relationships evolve from the diffuse ISM to molecular clouds.
An important parameter in this study is the total amount of dust along the line of sight, or dust column density, which is usually quantified by the dust optical depth τ (at 353 GHz). Because dust emission is optically thin at this frequency, this relates the modified blackbody (MBB) model of the emission to the total intensity via
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 dependent on the environment. This is not important for our results below and we use the MBB parameters τ and T_{d} from Planck Collaboration Int. XLVIII (2016), the calibration from Planck Collaboration XI (2014) at 353 GHz,
and the relation N_{H} = 5.87 × 10^{21} cm^{−2} × E(B − V)_{τ} (Bohlin et al. 1978; Rachford et al. 2009) to estimate N_{H}. It is preferable to use τ converted to N_{H} rather than an estimate of the gas column density derived from H I because of the presence of dust in the WIM that is sampled by all of our polarized and unpolarized observables.
We note that in this section we use not only the GNILC maps at 80′ and 160′ resolution, but also the alternative Stokes maps (Sect. 2.4) at finer resolutions of 40′, 20′, and 10′.
5.1. Origin of the observed variations of the polarization fraction p
The mutual correlations between p, 𝒮, and the column density N_{H} were studied in detail for the particular case of the Vela C molecular cloud by Fissel et al. (2016) using BLASTpol data. From the present Planck 2018 data, Fig. 12 shows how these correlations appear for the more diffuse ISM (4 × 10^{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. Colouring^{16} the data with 𝒮 (on a logarithmic scale) we see from the stratification of the data in Fig. 12 that there is a gradient of 𝒮 mainly perpendicular to the observed trend of p with respect to N_{H}. This analysis indicates that the decreases in p with 𝒮 and with N_{H} are mostly independent of each other.
Fig. 12. Polarization fraction p as a function of the total gas column density, N_{H}, coloured by the polarization angle dispersion function 𝒮 (on a logarithmic scale). The resolution of the data is 160′, in order to limit the bias in 𝒮. 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. 
Figure 13 shows how the variations of p and 𝒮 at a given column density are related to the dust temperature T_{d}. Dust tends to be systematically cooler when p is high and warmer when p is low (top panel). This is observed at all but the highest column densities (N_{H} > 8 × 10^{21} cm^{−2}). On the other hand, the opposite is seen in 𝒮 (bottom panel). The mirror symmetry between the two panels of Fig. 13 shows convincingly that there is in fact no physical relation between the polarization fraction p and the dust temperature T_{d} in the diffuse ISM. Even if it seems that, at any column density, high p corresponds to colder dust and low p to warmer dust, the bottom panel demonstrates that the value of p is actually driven by 𝒮, i.e., by the magnetic field structure and the depolarization produced by its variations along the line of sight and within the beam.
Fig. 13. Polarization fraction p (top panel), and polarization angle dispersion function 𝒮 (bottom panel) coloured by the dust temperature, T_{d}, as a function of the column density N_{H}. We note that the resolution of the data used here is 40′. Estimates of 𝒮 are nevertheless not biased, because only lines of sight with N_{H} ≥ 10^{21} cm^{−2} are considered. 
5.2. Exploring beyond firstorder trends using 𝒮 × p
In Sect. 4.3, we concluded that the inverse relationship between 𝒮 and p is a generic statistical property that results from the topology of the magnetic field alone, and that a trend 𝒮 ∝ 1/p, close to that observed, is expected on the basis of simple analytical models (Appendix F). It is therefore interesting to explore beyond this underlying cause for the inverse relationship, in search of evidence for the impact of other physical factors, such as dust alignment efficiency, elongation, and composition. For this we can use the product 𝒮 × p, which removes the impact of the magnetic field structure statistically. This does not mean that the product depends only on properties of the dust, e.g., the maximum polarization fraction p_{max}. As explained in Appendix F, the product 𝒮 × p also depends on the length over which dust structures are probed along the line of sight, as well as on the ratio of the turbulent to ordered magnetic field. Nevertheless, it is interesting to try this approach, as also emphasized by the mirror symmetry seen between the two panels of Fig. 13.
Accordingly, Fig. 14 compares the variations of not only p and 𝒮, but also 𝒮 × p with N_{H}, Galactic latitude b, and Galactic longitude l. It should be noted that throughout this entire analysis lines of sight close to the Galactic plane (b < 5°) are excluded.
Fig. 14. Twodimensional histograms with background colours encoding the density of points on a logarithmic scale, showing p (top), 𝒮 (middle), and 𝒮 × p (bottom) as a function of the column density N_{H} (left), Galactic latitude b (middle), and Galactic longitude l (right). The resolution is 160′. The colour bar shown in the top left panel is common to all plots. Black curves show the running means calculated as in Fig. 12, with error bars representing the scatter in each bin. For 𝒮, which is on a logarithmic scale, the median trend shown (thin black line) follows the density of points more faithfully than does the mean (thicker black line). 
As expected, the product 𝒮 × p has smaller variations with N_{H}, b, and l than exhibited by p and 𝒮 separately, and the decrease in 𝒮 × p with N_{H} is systematic, without significant departures.
Away from the Galactic plane, the dependence of 𝒮 × p on b is less pronounced than it is for p and also more symmetric between positive and negative latitudes. The strong dependence of p on b that can be attributed to the systematic change in the orientation of the mean magnetic field with respect to the line of sight is mitigated in 𝒮 × p, confirming our interpretation. However, there are still small variations of 𝒮 × p over a large spatial scale that remain to be interpreted. Towards the Galactic plane there is a pronounced dip, that is probably due to the accumulation of variously polarized structures along the line of sight at these low latitudes (Jones et al. 1992). This dependence on the latitude will be further discussed in Sect. 5.4.
As with the dependence on N_{H}, the variations of the product 𝒮 × p with Galactic longitude l are much less pronounced (of the order of 30%) than those of p and 𝒮 independently (a factor 3 or so in each case).
5.3. Dedicated study for six molecular regions in the Gould Belt
The radiative torques theory (RAT; Lazarian & Hoang 2007; Hoang & Lazarian 2016) makes strong predictions for the dependence of the dust alignment on local physical conditions, namely the intensity and anisotropy of the radiation field, and the angle between the magnetic field and the anisotropic radiation field. Dense regions, screened from the interstellar radiation field and possibly with embedded sources, should be promising regions in which to identify evidence for RATs (Vaillancourt & Andersson 2015; Wang et al. 2017).
To probe this possibility, we have selected a set of six 12° × 12° molecular regions in the Gould Belt (Dame et al. 2001). These are listed in Table 2. All but one (Aquila Rift) were already studied using Planck 2013 data in Planck Collaboration Int. XX (2015). The higher S/N in these bright, high column density regions enables an analysis at a higher resolution (40′, N_{side} = 256), and the uncertainty on the offset in total intensity I can be safely ignored. For this study, we therefore make use of the ASMs (Sect. 2.4) with the fiducial total intensity offset.
Selected molecular regions in the Gould Belt.
Figure 15 (top panel) shows that the variation of p with N_{H} is very diverse in these molecular clouds, as was already observed in Planck Collaboration Int. XX (2015). For some (e.g., Aquila Rift and Chamaeleon) p is fairly high in the more diffuse envelope but progressively decreases towards denser parts. Others (e.g., Polaris and Orion) have low p at all column densities, while for one (Ophiuchus) p increases, then decreases. In each region, the corresponding variations of 𝒮 with N_{H} (Fig. 15, middle panel) are clearly inversely related to those of p, so that the product is by contrast almost constant and uniform across the sample of clouds, as can be seen in Fig. 15 (bottom panel).
Fig. 15. Means from 2dimensional distributions of polarization properties and column density N_{H}, for selected regions in the Gould Belt, at a resolution of 40′: p (top); 𝒮 (middle); and 𝒮 × p (bottom). All bins in N_{H} contain the same number of pixels n, approximately 250. Error bars correspond to the uncertainty on the mean, i.e., , where σ is the statistical dispersion in the corresponding bin. The dashed horizontal line in the bottom panel is the mean value of 𝒮 × p at 160′ (Fig. 10), corrected for its dependence on the resolution, as per Eq. (10). 
Figure 16 (top panel) directly shows this inverse trend between 𝒮 and p in the selected molecular clouds. The various curves show the mean 𝒮 in each bin of p, all bins containing the same number of pixels. Also plotted is the average curve over all the different regions. Despite differences in mean column densities, regions as different as Polaris and Orion all fall on the same correlation line. This figure demonstrates that most of the variations of p with N_{H} in Fig. 15 can be attributed to variations of 𝒮 alone, and eventually to the variation of the magnetic field structure along the line of sight and in the plane of the sky.
Fig. 16. Mean 𝒮 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 𝒮 = f(p) trend at 160′ (Fig. 10), corrected for its dependence on the resolution, as per Eq. (10). All bins in p contain the same number of pixels, n ≈ 250. Error bars correspond to the uncertainty on the mean, i.e., , where σ is the statistical dispersion in the corresponding bin. 
We note that the inverse relationship between 𝒮 and p is the same as the one found in Sect. 4.2.2 over the full sky (Fig. 10), which is dominated by the highlatitude diffuse ISM, once the difference in resolution, and therefore in the lag δ used to compute 𝒮, is accounted for in the framework of our analytical model (see Appendix F.7). Indeed, for this analysis towards selected molecular regions, we work at a finer 40′ FWHM resolution, instead of 160′, and numerical results of the model show that the product 𝒮 × p scales as FWHM^{0.18} (Eq. (10)). The prediction of this model is shown as the dashed lines in Figs. 15 and 16.
The bottom panel of Fig. 16 shows the result of the same procedure applied to the phenomenological model described in Sect. 4.3 and smoothed to the same resolution of 40′. Our model, which was able to reproduce the trend 𝒮 ∝ 1/p observed at large scale (Figs. 10 and 11), can also reproduce it at the smaller scales probed here, in regions of 12° × 12°. The downward shift of the correlation observed in the data (compare Fig. 16 with Figs. 10 and 11), that is due to the change in resolution and is already integrated in our expression for 𝒮 × p, is also observed in the simulation.
As quantified in Appendix F within the framework of the phenomenological model of Planck Collaboration Int. XLIV (2016), the mean value of the 𝒮 × p product depends on f_{M}, the ratio of the dispersion of the turbulent component of the field to the mean field strength (see Eq. (F.53), where f_{m} scales linearly with f_{M}). Thus, the alignment of the data lines in the top panel of Fig. 16 is a remarkable result, which suggests that the strength of the turbulent component of the magnetic field, relative to the mean field strength, is comparable among Gould Belt clouds, and between clouds and the diffuse ISM, despite differences in the local starformation rate. This interpretation is illustrated by the correspondence between the top and bottom (data and model) plots of Fig. 16.
In the cold neutral medium, the magnetic and turbulent kinetic energies are known from H I Zeeman observations to be in approximate equipartition (Heiles & Troland 2005). Our analysis of the Planck data suggests that this energy equipartition also applies to the Gould Belt clouds. This result is consistent with the much earlier results derived from the modelling of stellar polarization data by Myers & Goodman (1991) and Jones et al. (1992).
To test the RAT theory, we need to estimate the relative intensity of the radiation field, G_{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 𝒮, i.e., look for a correlation between G_{0} and 𝒮 × p. To this end, we use an estimator G_{ℛ} of the radiation field intensity (Guillet et al. 2018; Fanciullo et al. 2015) that is based on the assumption of thermal equilibrium for large dust grains, which dominate the emission at this frequency. Under this hypothesis, the dust radiance ℛ, which is the integrated intensity of the dust emission (Planck Collaboration XI 2014), is balanced by the heating of dust by absorption of the ambient radiation field. The relative intensity G_{0} of this ambient field is therefore estimated using the radiance per unit visual extinction A_{V}, and in practice, the estimate G_{ℛ} is computed through
where E(B − V)_{ℛ} stands for the dust radiance converted to a colour excess using a correlation with extinction to quasars (Planck Collaboration XI 2014), while E(B − V)_{Green} is from a colour excess map (Green et al. 2018) based on PanSTARRS1 (PS1).
Figure 17 shows the correlation between this estimate G_{ℛ} of the radiation field intensity and the dust temperature T_{d} from the MBB fit, in the molecular regions selected, with data points coloured according to the dust optical depth τ converted to a column density N_{H}. This is the equivalent of Fig. 2 in Guillet et al. (2018), limited to the Gould Belt regions selected. We see that the correlation is quite tight and follows a scaling , 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, so that G_{ℛ} is overestimated. We also note that at high optical depths the spectral shape of the ambient radiation field is altered by the frequency dependence of the extinction, which would also impact the amount of energy absorbed and thus the dust temperature. Overall, over a wide range of G_{ℛ} this plot demonstrates that a change in the dust temperature is a good indicator of a change in the ambient radiation field intensity G_{0} in these molecular regions.
Fig. 17. Correlation between dust temperature T_{d} and our estimate G_{ℛ} 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). 
According to the RAT theory, we would expect grains in a more intense radiation field to be better aligned, and therefore that 𝒮 × p would tend to increase with increasing T_{d}. However, Fig. 18 does not show any correlation between the polarization fraction and the dust temperature, and the product 𝒮 × p, which we use as a proxy for dust alignment, does not show any trend with T_{d} either. This analysis for molecular clouds at a resolution of 40′ confirms our conclusion drawn in Sect. 5.1 for the diffuse ISM that there is no strong indication of a link between the polarization fraction p and the dust temperature T_{d}.
Fig. 18. Polarization fraction at 353 GHz (top) and product 𝒮 × 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 𝒮 × p at 160′ (Fig. 10), corrected for its dependence on the lag (Eq. (10)). 
5.4. Multiresolution view of the variations of p, 𝒮, and 𝒮 × p across the ISM
The discussion of the variations of p, 𝒮, and 𝒮 × p in the diffuse ISM at 160′ (Sect. 5.2) and in molecular clouds at 40′ (Sect. 5.3) suggests a multiresolution exploration of these trends. In Fig. 19, we present such a view by compiling mean trends at all resolutions, from 10′ to 160′. The impact of noise bias at low column densities in p, and even more so in 𝒮, is clearly seen as a rising deviation (dotted line) from the bundle of curves that otherwise roughly match – except for a global shift – at higher column densities, despite the fact that both p and 𝒮 have been debiased. All debiasing methods indeed fail when the S/N becomes lower than 1. This occurs below a threshold in N_{H} that is different for p and 𝒮 and increases with decreasing FWHM. Our debiasing methods are known to be statistically efficient when the S/N is higher than about 3 (Simmons & Stewart 1985; Montier et al. 2015a).
Fig. 19. Mean of p (top), 𝒮 (middle), and 𝒮 × 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 colour represents the density of points at a resolution of 10′, shown only for N_{H} > 4 × 10^{21} cm^{−2}. 
Figure 20 presents the evolution of the mean S/N in p and 𝒮 as a function of the column density for various resolutions of the map (still considering only b > 5°). Obviously, the S/N tends to increase with N_{H} at a given resolution, and to decrease with increasing resolution at a given N_{H}. With these figures we can estimate, at each resolution, the column density threshold above which the mean S/N for p and 𝒮 is larger than 3, i.e., above which debiased values of p and 𝒮 are robust. The S/N is smaller for 𝒮 than it is for p, and we note the large scatter in S/N for any given bin in N_{H}. Users of Planck data in polarization should take into account these thresholds in column density to estimate the reliability of the debiasing.
Fig. 20. Mean S/N of p (top), and 𝒮 (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, not to the uncertainty on the mean. The dashed line indicates the minimal S/N that ensures reliable mean values for debiased quantities. 
In Fig. 19, the data points below these thresholds are excluded from our analysis (dotted curves). The spread of values for p and 𝒮 × p at low column densities, induced by the uncertainty on the offset in total intensity I, is indicated by dashed lines for 160′. For reference, the density of points at a resolution of 10′ is plotted as a background. It is only plotted for N_{H} > 4 × 10^{21} cm^{−2}, which is the column density threshold for 𝒮 at this resolution.
Inspection of Fig. 19 shows that there always exists a range in column density where the curves at two consecutive resolutions are parallel to each other, i.e., they probe the same variations. It is therefore possible to obtain a unique, smooth and continuous trend for each quantity as a function of N_{H} through a renormalization of the profiles at each resolution, leading to Fig. 21. We proceed as follows, for p, 𝒮, and 𝒮 × p. At each resolution ω (expressed in arcmin), using Fig. 20 we determine the minimal column density above which 𝒮 is correctly debiased, and define a reference interval (indicated by the horizontal colour bars in Fig. 21). For two consecutive resolutions ω and ω/2, we compute the mean values of p, 𝒮, and 𝒮 × p at both resolutions on the common interval ℐ(ω)∩ℐ(ω/2) and then the ratio of these two values, r_{ω, ω/2}. Finally, we compute the factor by which each curve at resolution ω from Fig. 19 must be multiplied to be normalized to the curve at the coarsest resolution, ω_{max} = 160′, i.e., r(ω) = r_{ωmax, ωmax/2} × ... × r_{4ω, 2ω} × r_{2ω, ω}. This renormalization removes the depolarization induced in p by the smoothing of the data, as well as the change of the lag δ with the resolution, as far as 𝒮 in concerned.^{17}
Fig. 21. Mean 𝒮, p, and 𝒮 × p as a function of N_{H}, combining results from Planck maps at optimal resolutions for all lines of sight above b > 5° (solid curves). For clarity, 𝒮 has been raised vertically by a factor of 2. Upper and lower dashed red 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 on the mean, , 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 thin black dashed (b > 10°) and dotted (b > 2°) curves. Horizontal coloured bars indicate for each resolution ω the column density interval ℐ(ω) used in the renormalization procedure (see text). The green band highlights a 25% decrease in 𝒮 × p with column density up to 2 × 10^{22} cm^{−2}. 
The mirrored similarity of each detailed variation in the logarithmic representation of p and 𝒮 in Fig. 21 clearly emphasizes the inverse relationship between these two quantities. In our multiresolution normalized representation of the variations of p with the column density, the mean value of p decreases by a factor 3−4 from the lowest column densities at high latitudes and a resolution of 160′, to the highest column densities probed here at a resolution of 10′. This strong decrease is almost entirely mirrored as an increase in 𝒮 by the same factor, demonstrating again that the decrease in p with the column density is mainly a consequence of the depolarization by the structure of the magnetic field.
The residual trend in 𝒮 × p is small, about a 25% decrease with column density from 10^{20} cm^{−2} to 2 × 10^{22} cm^{−2}. For the case b > 10°, the profile of 𝒮 × p over this same range of N_{H} is quite flat. For the case including b > 2°, 𝒮 × p decreases systematically with N_{H}. In our phenomenological model, a systematic decrease in 𝒮 × p is expected at low Galactic latitudes from an increased number of independent layers N along the line of sight (see Eq. (F.45)), related to the increased dust opacity and/or length probed (Jones et al. 1992), that is not compensated by the inverse effect of increased 𝒮 due to an increased distance to the observed dust material (recall that 𝒮 depends on the physical scale probed in the cloud, therefore on its distance). There remains therefore little room for a systematic variation of grain alignment for column densities up to 2 × 10^{22} cm^{−2}.
At slightly higher column densities (N_{H} > 3 × 10^{22} cm^{−2}), we observe a decrease in p, together with an increase in 𝒮 and 𝒮 × p. Such a combination cannot be explained by a decrease in grain alignment, which would not affect 𝒮. These lines of sight, part of the Orion and Ophiuchus regions, are situated at intermediate latitudes (10° < b < 20°) and probably do not suffer from depolarization by the Galactic background, unlike other lines of sight at lower column densities situated at lower latitude. As can be seen from the density of points in Fig. 19, this departure from the mean trend has a low statistical significance, which prevents us from commenting further.
To conclude, most variations of the polarization fraction p with N_{H} are inversely related to those of 𝒮, a tracer of the depolarization by the turbulent magnetic field. The multiresolution study of the variation of 𝒮 × p with N_{H} does not indicate any systematic variation of the grain alignment efficiency beyond around 25%, up to a column density of 2 × 10^{22} cm^{−2}.
5.5. Grain alignment efficiency in the ISM
In this section, we discuss the impact of our results on the question of grain alignment.
Since the pioneering work of Myers & Goodman (1991) and Jones et al. (1992), it has been clear that the structure of the magnetic field along the line of sight plays a major role in the buildup of polarization observables. Nevertheless, the decrease in the polarization fraction with increasing column density is widely considered as evidence for a systematic decrease in the degree of grain alignment efficiency with increasing exinction (Whittet et al. 2008; Cashman & Clemens 2014; Alves et al. 2014).
In this paper, we have shown that most (if not all) variations observed in the polarization fraction p are mirrored in the dispersion of polarization angles, 𝒮, a quantity that is independent of the grain alignment efficiency and of dust optical properties. Quantitatively, the near constancy of 𝒮 × p with increasing column density indicates that the variations of the polarization fraction are dominated by the structure of the magnetic field, not only in the diffuse ISM, but also in molecular clouds, at least up to a column density of N_{H} ≈ 2 × 10^{22} cm^{−2}. The decrease in grain alignment efficiency with column density cannot exceed about 25%, from the most diffuse ISM up to this same column density, N_{H} ≈ 2 × 10^{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 fullsky data, we could not find, either in the diffuse ISM or in molecular clouds, any signature in polarization observables that could point to a significant variation of grain alignment related to a variation in the grain temperature.
However, the low resolution of Planck data (5′), combined with the smoothing of the maps necessary to guarantee that p and 𝒮 are not biased by noise (160′ in the high latitude diffuse ISM, 40′ in molecular clouds), does not allow us to probe the same physical conditions as, for example, NIR polarimetry through dense clouds (Jones et al. 2015). A detailed analysis of Planck polarization without the additional smoothing, and hence at higher resolution (7′) and higher column densities, will be pursued in a future paper dedicated to cold cores.
6. Comparison with starlight polarization at high Galactic latitudes
In this section, we correlate Planck polarization with starlight polarization in the diffuse ISM at high Galactic latitudes, and derive new constraints on dust models. Following the approach in Planck Collaboration Int. XXI (2015) for translucent lines of sight (0.15 < E(B − V)_{τ} < 0.8) at low Galactic latitudes (b < 30°), we can now derive emissiontoextinction polarization ratios for the diffuse highlatitude Galactic ISM (E(B − V)_{τ} < 0.15, corresponding to column densities N_{H} < 10^{21} cm^{−2}). The ratios are
Here P and I are what Planck has measured in the submillimetre. In the optical V band, p_{V} is the degree of polarization for a star to which the optical depth is τ_{V}. The latter is estimated from the colour excess of the star, E(B − V)^{⋆}, using A_{V} = R_{V} × E(B − V)^{⋆} with the ratio of total to selective extinction R_{V} = 3.1 (Fitzpatrick & Massa 2007), and τ_{V} = A_{V}/1.086.
These polarization ratios quantify the amount of polarized emission per unit polarized extinction. Because they measure the effects of the same grains at different wavelengths, they remove the firstorder dependencies of the polarization observables on the magnetic field structure and grain alignment efficiency (Planck Collaboration Int. XXI 2015). As such, they directly provide observational constraints on the optical properties of the dust population that is aligned with the magnetic field, and strongly constrain dust models (Guillet et al. 2018).
The second of these ratios, R_{S/V}, being inversely proportional to I, is sensitive to the total intensity offset. We comment on the derived values of R_{S/V} for the low, fiducial, and high values of the offset in Sects. 6.4 and 6.5.
We are interested in examining how the amount of submillimetre polarization is related to the amount of optical polarization from the same dust. If the dust probed in the submillimetre and the optical is the same, then the polarization orientations should be orthogonal. As discussed in Appendix G.4, this is the case for the lines of sight used in our analysis of the ratios. This is quantified by the polarization angle difference Δψ_{S/V} which takes into account the 90° difference (see Eq. (G.6)).
In Appendix G we describe in detail our estimates of the noise and systematic uncertainties that affect the submillimetre and optical data as used in this new analysis. We highlight the relevant results in the discussion below.
6.1. Estimates for starlight reddening
To enable appropriate comparison of polarization properties in the optical and in the submillimetre on lines of sight to stars, it is necessary to obtain an estimate of the reddening to the star. To this end, we obtain the distance D^{⋆} to each star by extracting the parallax θ^{⋆} from the Gaia DR2 release (Gaia Collaboration 2016, 2018) or from the polarization catalogues when Gaia data are not available. Then we derive an estimate of the reddening to the star, E(B − V)^{⋆}, by interpolating the PS1based 3dimensional reddening data cube (Green et al. 2018) at distance D^{⋆}. This data cube is composed of 31 maps, each representing a range in distance modulus. To limit the impact of noise in our analysis, the 31 maps were separately smoothed to a resolution of 40′ and downgraded to N_{side} = 256. We also converted the PS1based reddening to the Johnson E(B − V) scale using the relation E(B − V)_{Johnson} = E(B − V)_{PS1}/1.0735 (Table 6 of Schlafly & Finkbeiner 2011). The Johnson scale is used hereafter without explicit subscripting. From the same maps we also obtain the total reddening along the line of sight, E(B − V)^{∞}. Uncertainties related to the reddening maps are estimated in Appendix G.3.
6.2. Polarization data
For this analysis, we use the alternative Planck 353GHz 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, 2004, 2014; Berdyugin & Teerikorpi 2001, 2002), we extract data for 2461 lines of sight to stars with measured degree of polarization p_{V}, uncertainty σ_{pV}, and polarization angle ψ_{V} (in the IAU convention, consistently with our definition of ψ for the Planck data). These catalogues cover the northern Galactic hemisphere at high latitudes (2092 stars with b > 30°), and part of the southern cap (369 stars with b < −59°). After removing 70 stars falling outside the region covered by PS1 (mainly in the southern hemisphere) and 3 stars without a distance estimate, there remain 2388 stars for which we have both reddening estimates and optical polarization data.
As with the Planck submillimetre data, we use the MAS estimator (Plaszczynski et al. 2014) to debias the degree of polarization p_{V} in the optical. Using these values of p_{V} and ψ_{V}, we then define Stokes parameters in extinction, q_{V} and u_{V}, in the same HEALPix convention as our Planck data:
Emissiontoextinction ratios are subject to systematic errors because extinction probes the ISM in the foreground to the star, while emission probes the entire line of sight (Planck Collaboration Int. XXI 2015). Figure 22 presents the histogram of the ratio E(B − V)^{⋆}/E(B − V)^{∞} of the reddening to the star to the total reddening, i.e., the fraction of ISM material that is in front of each star. The median ratio for our full sample is 0.83, illustrating the potential for systematic effects on the polarization ratios.
Fig. 22. Histogram of the ratio of the reddening to the star to the total reddening on the same line of sight, E(B − V)^{⋆}/E(B − V)^{∞}, as derived from the PanSTARRS1 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. 
If we assumed for simplicity that the ISM along the line of sight were uniform (in density, magneticfield orientation, and dust properties), then the polarization ratio R_{P/p} = P/p_{V} would artificially increase linearly with decreasing E(B − V)^{⋆}/E(B − V)^{∞}. Consequently, by neglecting the presence of background material we would typically overestimate the polarization ratio R_{P/p} by 17%. Given this contamination, to debias our estimate of R_{P/p} we replace p_{V} by a linear estimate of what its value would be if the star were at infinity:
with an associated uncertainty
We use similar expressions to estimate and and their uncertainties and . On the other hand, R_{S/V}, as a ratio of nondimensional quantities, would be unaffected by this uniform background.
6.3. Selection of the lines of sight
However, the ISM is not uniform and as a consequence our estimates of both R_{S/V} and R_{P/p} could be biased by the presence of a background whose properties are different from those of the foreground to the star. The magnitude of this uncertainty is evaluated in Appendix G.2. We minimize the contribution of this uncertainty by excluding those lines of sight with an important background, as inferred from the ratio E(B − V)^{⋆}/E(B − V)^{∞} shown in Fig. 22. We explicitly choose to keep only stars for which E(B − V)^{⋆}/E(B − V)^{∞} > 0.75.
We also exclude those lines of sight where Δψ_{S/V} is significantly different than the expected 0°, having found that for these lines of sight the rms scatter about the best fit correlations yielding R_{P/p} and R_{S/V} (Sect. 6.4) is indeed much larger. To be conservative and retain enough lines of sight for our subsequent statistical analysis, we excluded only lines of sight for which Δψ_{S/V} > 45°, i.e., about 9% of the sample. Our results are not sensitive to this particular choice.
Our final sample contains 1505 stars. The ISM towards these stars in emission is representative of diffuse dust at high Galactic latitudes, with MBB fit parameters T_{d} = 19.7 ± 1.3 K, β = 1.60 ± 0.15, and E(B − V)_{τ} ∈ [0.01, 0.24], with a mean reddening ⟨E(B − V)_{τ}⟩ = 0.03.
6.4. Determination of the polarization ratios
Following the approach in Planck Collaboration Int. XXI (2015), we derive R_{P/p} through a joint correlation of the pair (Q, U) with (), and derive R_{S/V} through a joint correlation of (Q/I, U/I) with (q_{V}/τ_{V}, u_{V}/τ_{V}). In Fig. 23 we present the two correlation scatterplots, that for R_{P/p} on the left, and that for R_{S/V} on the right. For the fitting, we settle on the Bayesian method of Kelly (2007), but we obtain the same results with other fitting methods making use of uncertainties on both axes. Both for determining the value of the ratio (the slope) and for calculating the reduced χ^{2} to assess the quality of the fit, it is important to assess all sources of uncertainty, as addressed in this section and in more detail in Appendix G.
Fig. 23. Correlation between Stokes polarization parameters in emission at 353 GHz and in optical extinction, with the colour in the 2dimensional histogram representing the density of points. Left: Stokes parameters (Q, U) versus (), 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). The reduced χ^{2}, the Pearson correlation coefficient, and the correlation coefficient inferred from the Bayesian method (Kelly 2007) are listed. 
As a first test, we fit the data for the 206 translucent lines of sight from Planck Collaboration Int. XXI (2015) with this estimator and find no change, even though we smooth Planck Stokes parameter maps to 40′ FWHM. We are therefore confident that it is legitimate to compare the polarization ratios that we derive here at 40′ resolution with those measured at 7′ resolution in Planck Collaboration Int. XXI (2015).
The correlation of (Q, U) with () shown in Fig. 23 (left) is tight, with a Pearson correlation coefficient −0.92. For lines of sight where p is low, error bars have been greatly increased by the correction factor for beam depolarization (Appendix G.1). When systematic uncertainties are taken into account, the reduced χ^{2} is good (1.02, compared to 3.9 when they are ignored). The fit yields a polarization ratio R_{P/p} = (5.42 ± 0.05) MJy sr^{−1}, similar to the value found for translucent lines of sight, (5.4 ± 0.5) MJy sr^{−1} (Planck Collaboration Int. XXI 2015).
We find a somewhat larger scatter in the correlation plot (Q/I, U/I) with (q_{V}/τ_{V}, u_{V}/τ_{V}) in Fig. 23 (right), as quantified by the slightly lower absolute value of the correlation coefficient. Here also, the reduced χ^{2} is good when systematic uncertainties are included (0.96), and larger when they are not (3.0). The fitted slope R_{S/V} = 4.31 ± 0.04 is also compatible with the value 4.2 ± 0.5 found for translucent lines of sight at 7′ resolution (Planck Collaboration Int. XXI 2015). We find similar values at other resolutions: R_{S/V} = 4.2 at 80′ resolution; and R_{S/V} = 4.2 at 20′ resolution. These results are for the fiducial intensity offset. With the low and high offsets at 40′ resolution, we obtain R_{S/V} = 4.8 and R_{S/V} = 3.9, respectively, which makes the Planck intensity offset the main source of uncertainty on R_{S/V}.
6.5. Variations of R_{P/p} and R_{S/V} with column density
Our sample contains enough lines of sight to study variations of the polarization ratios with column density N_{H}, which is determined from the dust optical depth at 353 GHz as explained at the beginning of Sect. 5. This is presented in Fig. 24. The sample is divided into 15 independent bins in N_{H}, each bin containing 100 stars. For the polarization ratio R_{P/p} at low N_{H}, we observe a roughly 10% increase with increasing N_{H}, from about 5.0 MJy sr^{−1} at N_{H} ≈ 6 × 10^{19} cm^{−2} to 5.4 MJy sr^{−1} at N_{H} ≈ 1.5 × 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. Both the Pearson correlation coefficients and the correlation coefficients provided by the Bayesian method (Kelly 2007) are high enough at all column densities to bring confidence in our results. For comparison, the values of the polarization ratios found for translucent lines of sight in Planck Collaboration Int. XXI (2015), together with their ranges of uncertainty, are also displayed in Fig. 24. Altogether, R_{P/p} and R_{S/V} are remarkably constant with column density and consistent with the values determined for translucent lines of sight. We note that there are some small variations of potential significance, such as the 10% increase in R_{P/p} over the range 6 × 10^{19} − 1.5 × 10^{20} cm^{−2}, but it is beyond the scope of this paper to discuss them.
Fig. 24. Emissiontoextinction 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}. For the running mean each bin contains the same number (100) of lines of sight. The lower dark blue dotteddashed lines indicate the reduced χ^{2} of the fits (with the scale on the left axis). Dashed red and purple curves represent the Pearson and Bayesian (Kelly 2007) 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). 
6.6. Maximum p_{V}/E(B − V) at low column densities
Regarding the observed starlight polarization fraction per unit reddening in the optical, p_{V}/E(B − V), its maximum value was first estimated to be at most 9% by Serkowski et al. (1975). This maximum value is often considered as providing an upper limit on the dust alignment efficiency, although it is based on less than 300 stars at moderate extinction (E(B − V) > 0.15), characteristic of translucent lines of sight. Several attempts have been made to constrain this maximum value at low reddening (Fosalba et al. 2002; Frisch et al. 2015; Skalidis et al. 2018), suggesting larger values, but with poor precision. Our present sample probes more diffuse lines of sight and, with more statistical significance, allows us to characterize the maximum polarization fraction at low reddening, extending to E(B − V) ≈ 0.01.
To study the dependence of the polarization fractions P/I and p_{V}/E(B − V) on column density, we combine data for the 1505 diffuse lines of sight to stars from this study (high latitude, low E(B − V), 40′ resolution for Planck data and E(B − V) map) and the 206 stars on translucent lines of sight (low latitude, moderate E(B − V), 7′ resolution for Planck data) from Planck Collaboration Int. XXI (2015).
Figure 25 shows how these polarization fractions in emission and extinction vary with column density. Polarization fractions at 353 GHz never reach low values because Planck polarization data have been corrected for beam depolarization through Eq. (F.63). However, this does not affect our results for the high percentiles because high values of p suffer very little depolarization. We overplot the upper limit p_{V} ≤ 9%E(B − V) of Serkowski et al. (1975), the nonlinear fit p_{V}/E(B − V) ∝ E(B − V)^{−0.2} proposed by Fosalba et al. (2002) for polarization in extinction (Fig. 25, right), and the estimate for the maximum value of the polarization fraction in emission observed for translucent lines of sight on the left (approximately 14%, Planck Collaboration Int. XXI 2015).
Fig. 25. Polarization fraction in emission at 353 GHz, p = P/I (left), and in optical extinction, p_{V}/E(B − V) (right), as a function of the column density, N_{H}. The sample in blue shows 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 (1505 stars), where data have been corrected for systematic effects such as beam depolarization (see Appendix G). 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) ∝ E(B − V)^{−0.2}, is shown for comparison. 
Assuming a maximum polarization in emission of p_{max} = 20% for our sample at a resolution of 40′ (close to its 99th percentile^{18}), and a polarization ratio R_{S/V} = 4.3 (Fig. 23), we would expect a maximum polarization fraction in extinction of p_{V}/E(B − V)_{max} = (3.1 × p_{max})/(R_{S/V} × 1.086)≈13%. As seen in Fig. 25 (right), this upper limit is somewhat smaller than the measured 99th percentile of our data in extinction. We would reach a similar conclusion using the value p_{max} = 22% that we obtained for the diffuse ISM in Sect. 4.1.1, with p_{V}/E(B − V)_{max} ≃ 14.5%. This upper limit is also smaller than the upper limit proposed by Fosalba et al. (2002) based on a study of the dependence of the mean starlight polarization fraction p_{V}/E(B − V) with E(B − V).
However, the distribution in the density of points for polarization in extinction (right panel) compared to that for polarization in emission (left panel) suggests that lines of sight with high p_{V}/E(B − V) might be outliers. One should indeed be aware that the direct derivation of the maximum polarization fraction in extinction is much more subject to noise and systematics than our derivation, which is based on the measurement of the polarization ratio R_{S/V} and the much better characterized maximum polarization fraction p in emission. We therefore consider the value of 13% as a wellconstrained 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%. Such a decrease is usually attributed to a loss of grain alignment (see Andersson et al. 2015 and references therein). However, inspection of Fig. 21 for Galactic latitudes higher than 10° shows that the product 𝒮 × p remains constant over the range of column densities probed here. Following our analysis in Sect. 5, we therefore attribute most of this decrease in the maximum polarization fraction when the column density increases to the increasing depolarizing effect from the structure of the magnetic field along the line of sight, with little room for a systematic decrease in the grain alignment efficiency.
Dust models should therefore be able to reproduce the maximum observed polarization fractions, p_{max} = 20% in emission and (p_{V}/E(B − V))_{max} = 13% in extinction, even when applied to regimes in column densities where such values are actually never directly observed.
7. Conclusions
In this paper, we have analysed the Planck 2018 thermal dust polarization data at 353 GHz. Starting from fullsky maps of Stokes I, Q, and U at a uniform 80′ resolution, processed with the Generalized Needlet Internal Linear Combination (GNILC) algorithm (Remazeilles et al. 2011) to mitigate the contamination by CIB and CMB anisotropies as well as noise, we have presented the maps of polarization fraction p, polarization angle ψ, and polarization angle dispersion function 𝒮, with their associated uncertainties. The statistical analysis of these maps provides new insights into the astrophysics of dust polarization.
We have been able to determine the maximum observed polarization fraction, , at this resolution and frequency, where the second uncertainty is statistical, underscoring the excellent quality of the data, and the first reflects the principal systematic uncertainty affecting this determination, which is linked to the uncertainty on the Galactic emission zero level in total intensity (Planck Collaboration XI 2014). This maximum polarization fraction provides strong constraints for models of dust properties and alignment in the Galactic magnetic field (Guillet et al. 2018). Owing to the strong effect of the magnetic field morphology, a low value for the maximum polarization fraction in a given region is not an indication that grain alignment in that region is ineffective, but rather that polarization is strongly affected by depolarization because the direction of the largescale field is closer to the line of sight in that region.
We confirmed that the statistical properties of p, ψ, and 𝒮 essentially reflect the structure of the Galactic magnetic field (Planck Collaboration Int. XIX 2015), with a clear peak of the polarization angle near 0°, corresponding to the field being parallel to the Galactic plane, and an inverse proportionality between the polarization fraction p and the polarization angle dispersion function 𝒮. We showed analytically, and using numerical models of the polarized sky, that this relationship can be reproduced statistically to first order by an interstellar magnetic field incorporating a turbulent component superposed on a small number of layers with a simple uniform field configuration. Looking for evidence in the diffuse ISM (N_{H} < 8 × 10^{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 𝒮, which does not depend on the dust physics.
Based on this analysis, we introduced the product 𝒮 × p as a means of exploring the nongeometric elements of the polarization maps, such as variations in grain properties, in alignment physics, or in ISM phase contributions. We showed that 𝒮 × p exhibits smaller and smoother variations than either p or 𝒮 when considered as a function of the Galactic latitude b, the Galactic longitude l, or the column density (which is simply scaled from the dust optical depth τ at 353 GHz).
We provided an analysis at a finer angular resolution of 40′ using the Planck 2018 data, towards a sample of six molecular regions in the Gould Belt. This confirmed the trends observed at coarser resolution over the full sky, most notably that the polarization angle dispersion function is inversely proportional to the polarization fraction, 𝒮 ∝ 1/p. Strikingly, the 𝒮 versus p curves for the different regions all fall on the same line, demonstrating that most of the variations of p with column density are driven by variations of 𝒮, i.e., by the structure of the magnetic field along the line of sight (Planck Collaboration Int. XX 2015). Considering then the product 𝒮 × p and how it varies with dust temperature T_{d} in these regions, we found no evidence for a link between the polarization properties and the intensity of the radiation field. Comparing these properties with column density in a multiresolution view, we found that the product 𝒮 × 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 of 3−4 over the same interval.
We also compared the Planck polarization data with optical stellar polarization data in the high Galactic latitude sky, expanding on the analysis done in Planck Collaboration Int. XXI (2015) for translucent lines of sight. We constrained the polarization properties of the dust at low column densities, quantified by the polarization ratios R_{P/p} = P/p_{V} and R_{S/V} = (P/I)/(p/τ_{V}) defined in Planck Collaboration Int. XXI (2015). The larger statistical sample (1505 stars selected) allowed us to study the dependence of these polarization ratios on column density. We found R_{P/p} to increase from 5.0 MJy sr^{−1} at N_{H} = 6 × 10^{19} cm^{−2} to 5.4 MJy sr^{−1} for N_{H} > 1.5 × 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.3) with that for translucent lines of sight (4.2 ± 0.5), having a decreasing, flat, or slightly increasing trend with the column density, depending on the offset for the Planck intensity map at 353 GHz, which is here again the dominant systematic effect. Combining emission and extinction measurements, we derived an estimate for the maximum value of the polarization fraction in the visible at high Galactic latitude, p_{V}/E(B − V) ≤ 13%, significantly higher than the value of 9% characterizing translucent lines of sight at low latitudes (Serkowski et al. 1975). This is a strong new constraint that dust models must satisfy.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
In this paper, all maps are shown either with a Mollweide projection of the full sky, in Galactic coordinates centred on the Galactic centre (GC), or with an orthographic projection of both hemispheres. In this latter case, the northern Galactic hemisphere is always on the left and the southern Galactic hemisphere on the right, with the rotation of each hemisphere such that the Galactic centre (l = 0°) is towards the top.
We note a typo in Eq. (B.3) of Planck Collaboration Int. XIX (2015), which should include a factor of 1/P on the righthand side.
When considering the Monte Carlo simulations discussed in the previous subsection, we find that the ratio of the ensemble average map ⟨𝒮⟩ to the map 𝒮 computed from the smoothed GNILC data have 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.
Other advanced diagnostics from polarization gradients are discussed in Herron et al. (2018), but further discussion of these is beyond the scope of this paper.
Those are shown as dashed curves in Fig. 5.
We note that the linear fitting of the mean log 𝒮 per bin of log p that was originally used in Planck Collaboration Int. XIX (2015) and Planck Collaboration Int. XX (2015) is not the optimal procedure to quantify the inverse relationship between 𝒮 and p.
This parameter is related to p_{max}, the maximum polarization fraction observed, by p_{max} = p_{0}/(1 − p_{0}/3) (Planck Collaboration Int. XX 2015).
We write E(B − V)_{τ} instead of simply E(B − V) to emphasize that this colour excess is computed from the dust optical depth derived from Planck data, and to distinguish it from other estimates used in Sect. 6.
For 𝒮 × p, this renormalization is consistent with the scaling with the resolution, 𝒮 × p ∝ ω^{0.18} (Eq. (10)).
The percentile curves in Fig. 25 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.
Because the polarization angle dispersion function 𝒮 involves an average over lags in [δ/2, 3δ/2], we note that f_{m}(δ) actually stands for an average of the fluctuation ratio of the magnetic field over this range (see Appendix F.7.)
Whether we add the Gaussian random realizations to the submillimetre or to the optical Stokes parameters does not affect the resulting histogram of differences in polarization angles and so we do not show the conversion of the equations for the submillimetre beam depolarization into their analogous form in the optical.
Acknowledgments
The Planck Collaboration acknowledges the support of ESA; CNES, and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planckcollaboration. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 Research & Innovation Framework Programme/ERC grant agreement ERC2016ADG742719. This research has received funding from the Agence Nationale de la Recherche (ANR17CE310022). We thank Pekka Teerikorpi and Andrei Berdyugin for kindly providing stellar polarization data and providing insights on stellar polarization references, Gina Panopoulou for discussions of the zero point calibration of the polarization angle, and Ralf Siebenmorgen and Nikolai Voshchinnikov for statistical discussions. We gratefully acknowledge the help of Rosine Lallement regarding the handling of the Gaia data.
References
 Alina, D., Montier, L., Ristorcelli, I., et al. 2016, A&A, 595, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alves, F. O., Frau, P., Girart, J. M., et al. 2014, A&A, 569, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alves, M. I. R., Boulanger, F., Ferrière, K., & Montier, L. 2018, A&A, 611, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andersson, B.G., & Potter, S. B. 2010, ApJ, 720, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Andersson, B.G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Aumont, J., MacíasPérez, J. F., Ritacco, A., Ponthieu, N., & Mangilli, A. 2018, A&A, 634, A100 [Google Scholar]
 Benoît, A., Ade, P., Amblard, A., et al. 2004, A&A, 424, 571 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugin, A., & Teerikorpi, P. 2001, A&A, 368, 635 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugin, A., & Teerikorpi, P. 2002, A&A, 384, 1050 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugin, A., Teerikorpi, P., Haikala, L., et al. 2001, A&A, 372, 276 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugin, A., Piirola, V., & Teerikorpi, P. 2004, A&A, 424, 873 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugin, A., Piirola, V., & Teerikorpi, P. 2014, A&A, 561, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhart, B., Lazarian, A., & Gaensler, B. M. 2012, ApJ, 749, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Cardoso, J., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, in IEEE Journal of Selected Topics in Signal Processing, special issue on Signal Processing for Astronomical and Space Research Applications, 2, 735 [Google Scholar]
 Cashman, L. R., & Clemens, D. P. 2014, ApJ, 793, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, L., & Greenstein, J. L. 1949, Phys. Rev., 75, 1605 [CrossRef] [Google Scholar]
 Davis, Jr., L., & Greenstein, J. L. 1951, ApJ, 114, 206 [Google Scholar]
 de Bernardis, P., Ade, P. A. R., Artusa, R., et al. 1999, Nature, 43, 289 [Google Scholar]
 Fanciullo, L., Guillet, V., Aniano, G., et al. 2015, A&A, 580, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2016, ApJ, 824, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fosalba, P., Lazarian, A., Prunet, S., & Tauber, J. A. 2002, ApJ, 564, 762 [NASA ADS] [CrossRef] [Google Scholar]
 Frisch, P. C., Berdyugin, A., Piirola, V., et al. 2015, ApJ, 814, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, PASA, 25, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Gaensler, B. M., Haverkorn, M., Burkhart, B., et al. 2011, Nature, 478, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration 2018, VizieR Online Data Catalog: I/345 [Google Scholar]
 Gautier, III., T. N., Boulanger, F., Perault, M., & Puget, J. L. 1992, AJ, 103, 1313 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, T., Boulanger, F., Martin, P. G., et al. 2017, A&A, 601, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gold, B., Odegard, N., Weiland, J. L., et al. 2011, ApJ, 192, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hall, J. S. 1949, Science, 109, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Hazewinkel, M. 2013, Encyclopaedia of Mathematics, Encyclopaedia of Mathematics (Netherlands: Springer) [Google Scholar]
 Heiles, C. 2000, AJ, 119, 923 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 [NASA ADS] [CrossRef] [Google Scholar]
 Herron, C. A., Gaensler, B. M., Lewis, G. F., & McClureGriffiths, N. M. 2018, ApJ, 853, 9 [CrossRef] [Google Scholar]
 Hildebrand, R. H. 1988, Astrophys. Lett. Commun., 26, 263 [Google Scholar]
 Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Hiltner, W. A. 1949, Science, 109, 165 [Google Scholar]
 Hoang, T., & Lazarian, A. 2016, ApJ, 831, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, T. J., Klebe, D., & Dickey, J. M. 1992, ApJ, 389, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, T. J., Bagley, M., Krejny, M., Andersson, B.G., & Bastien, P. 2015, AJ, 149, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368 [NASA ADS] [CrossRef] [Google Scholar]
 Kelly, B. C. 2007, ApJ, 665, 1489 [NASA ADS] [CrossRef] [Google Scholar]
 Kruegel, E. 2003, The Physics of Interstellar Dust (Bristol, UK: The Institute of Physics) [CrossRef] [Google Scholar]
 Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910 [Google Scholar]
 Levrier, F., Neveu, J., Falgarone, E., et al. 2018, A&A, 614, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montier, L., Plaszczynski, S., Levrier, F., et al. 2015a, A&A, 574, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montier, L., Plaszczynski, S., Levrier, F., et al. 2015b, A&A, 574, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Myers, P. C., & Goodman, A. A. 1991, ApJ, 373, 509 [NASA ADS] [CrossRef] [Google Scholar]
 NaghizadehKhouei, J., & Clarke, D. 1993, A&A, 274, 968 [NASA ADS] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJ, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration 2018, The Legacy Explanatory Supplement (ESA), http://wiki.cosmos.esa.int/plancklegacyarchive [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2020, A&A, 641, A2 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2020, A&A, 641, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2020, A&A, 641, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2020, A&A, 641, A9 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2020, A&A, 641, A10 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2020, A&A, 641, A12 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XVII. 2014, A&A, 566, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXI. 2015, A&A, 576, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXIX. 2016, A&A, 586, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXII 2016, A&A, 586, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXIII. 2016, A&A, 586, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXIV. 2016, A&A, 586, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXV. 2016, A&A, 586, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLII. 2016, A&A, 596, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLIV. 2016, A&A, 596, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVI. 2016, A&A, 596, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plaszczynski, S., Montier, L., Levrier, F., & Tristram, M. 2014, MNRAS, 439, 4048 [NASA ADS] [CrossRef] [Google Scholar]
 Rachford, B. L., Snow, T. P., Destree, J. D., et al. 2009, ApJ, 180, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011, MNRAS, 418, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Savage, A. H., Spangler, S. R., & Fischer, P. D. 2013, ApJ, 765, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, ApJ, 196, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Simmons, J. F. L., & Stewart, B. G. 1985, A&A, 142, 100 [NASA ADS] [Google Scholar]
 Skalidis, R., Panopoulou, G. V., Tassis, K., et al. 2018, A&A, 616, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, J. D., Bracco, A., & Pon, A. 2018, A&A, 609, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, W. 1966, ApJ, 144, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90 [NASA ADS] [CrossRef] [Google Scholar]
 Vaillancourt, J. E., & Andersson, B.G. 2015, ApJ, 812, L7 [CrossRef] [Google Scholar]
 Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, J.W., Lai, S.P., Eswaraiah, C., et al. 2017, ApJ, 849, 157 [CrossRef] [Google Scholar]
 Whittet, D. C. B., Hough, J. H., Lazarian, A., & Hoang, T. 2008, ApJ, 674, 304 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [Google Scholar]
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 353GHz polarized sky at 1° resolution, focusing on the statistics of the polarization fraction p and polarization angle ψ, at low and intermediate Galactic latitudes. We found a high maximum polarization fraction (p_{max} > 19.8%) in the most diffuse regions probed. This maximum polarization fraction was found to decrease as total gas column density N_{H} increases, which might be interpreted as changes of grain alignment properties or as the effect of magneticfield structure along the line of sight. We also characterized the structure of the polarization angle map by computing its local dispersion function 𝒮, which was found to be inversely related to the polarization fraction, lending support to the second explanation.
In Planck Collaboration Int. XX (2015), we presented an analysis of Planck 353GHz 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 these constraints into account.
In Planck Collaboration Int. XXXII (2016), we studied the correlation between the magneticfield 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 magneticfield geometry of interstellar filaments in Planck 353GHz 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 in polarization fraction within the filaments, although this could be due not only to the effect of magnetic field tangling along the line of sight, but also in part to changes of grain alignment properties.
In Planck Collaboration Int. XXXIV (2016), we combined the polarization data from Planck at 353 GHz with rotation measure (RM) observations from Savage et al. (2013) towards a massive starforming region, the Rosette Nebula in the Monoceros molecular cloud, to study the impact of an expanding H II region on the morphology of the Galactic magnetic field. We introduced an analytical model that describes the magnetic field structure in a spherical shell, following the expansion of an ionized nebula in a medium with uniform density and magnetic field, and fitted it to the data. This work was subsequently extended to nonspherical bubbles to model the structure of the largescale magnetic field in the Local Bubble (Alves et al. 2018). The Planck polarization data towards the Orion–Eridanus superbubble provide additional evidence for the impact of massive stars on the magnetic field structure (Soler et al. 2018).
In Planck Collaboration Int. XXXV (2016), we studied the relative orientation between filamentary structures of matter and the magnetic field towards molecular clouds of the Gould Belt, probing lines of sight with total gas column densities N_{H} from around 10^{21} 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 N_{H} of a few 10^{21} cm^{−2}. This change in relative orientation was found to be compatible with simulations of sub or transAlfvé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 Hessianextracted filaments at angular scales where the E/B powerasymmetry and TE correlation are observed (Planck Collaboration Int. XXX 2016). The preferred orientation of these filaments parallel to the magnetic field is able to account for both of these observations and was quantified by an oriented stacking of the maps of I, Q, U, E, and B. From these stacked maps and the histogram of relative orientations, we derived an estimate of the mean polarization fraction in the filaments, ⟨p⟩≈11%.
In Planck Collaboration Int. XLII (2016), we provided a comparison of three different models of the largescale 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 using the 1point statistics of p and ψ, we could constrain the orientation of the largescale 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 of the strengths of the turbulent to mean components of the GMF (f_{M} = 0.9 ± 0.1). This phenomenological framework was further improved by Ghosh et al. (2017) and Vansyngel et al. (2017).
Appendix B: GNILC Stokes and covariance maps
B.1. Variable resolution GNILC maps
For reference, in Fig. B.1 we show the GNILC Stokes maps at 353 GHz and variable resolution over the sky, alongside the map of the effective FWHM, whose discrete values are 5′, 7′, 10′, 15′, 20′, 30′, 60′, and 80′. The total intensity map corresponds to the fiducial offset value.
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′. The scale of the I map is logarithmic, while the rest are linear. 
B.2. GNILCprocessed 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 mimimumvariance 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^{19}. The residual noise rms fluctuations, N_{GNILC}, of the GNILC dust map are thus related to the instrumental noise rms fluctuations, N_{i}, in each Planck frequency map X_{i}, according to
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 halfmission maps X_{i, HM1} and X_{i, HM2}:
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 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 multiresolution case. The estimate 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}^{20}, 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 times the resolution of the Stokes maps. We then compute the local (co)variance value in each region of a given resolution of the Planck and GNILC noise maps, and , 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 scaled according to
The resulting covariance maps σ_{XY, GNILC} (p) are what we refer to as the GNILCprocessed covariance maps in the rest of this paper.
We note that they are built using the PSBonly data for polarization at 353 GHz, both for the uniform 80′ resolution case and for the Bmodedriven, 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 to MJy^{2} sr^{−2} using the conversion factor at 353 GHz (Planck Collaboration III 2020).
Figure B.2 shows these covariance maps for the variable resolution case, while Fig. B.3 shows the covariance maps at the common, uniform 80′ resolution. The sky patterns of these uniform resolution covariance maps are by construction extremely similar to those directly taken from the Planck 2018 data, but with a significant improvement in the amplitudes, as shown in Table B.1, which gives the mean ratios σ_{XY, GNILC}/σ_{XY} of the GNILC covariance maps σ_{XY,GNILC} at 80′ resolution and N_{side} = 128 to the corresponding maps σ_{XY} in the Planck 2018 data release. Consequently, we use these GNILCprocessed covariance maps to assess statistical uncertainties in our analysis.
Fig. B.2. GNILCprocessed covariance maps at 353 GHz and varying resolution. From top to bottom and left to right, they are σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU}. 
Fig. B.3. GNILCprocessed covariance maps at 353 GHz and unifom 80′ resolution. From top to bottom and left to right, they are σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU}. 
GNILC versus Planck 2018 data release covariances.
Appendix C: Endtoend simulations
The quality of the data presented here is assessed through a series of endtoend (E2E) simulations that take into account all the known systematics and noise properties of the data. The process begins with a model of the sky (including foregrounds, in both total intensity and polarization), from which timelines (including all known effects) are simulated. These timelines are then processed through the Planck mapmaking pipeline (Planck Collaboration III 2020). These E2E simulations are the same ones used in Planck Collaboration XI (2020) and we refer the reader to appendix A of that paper for a detailed description. The dust component of the model used in these simulations is a combination of a realization of the Vansyngel et al. (2017) statistical model for ℓ ≥ 20 and the actual Planck 2018 Q and U maps at 353 GHz for ℓ ≤ 10, with a smooth transition between the two in the 10 ≤ ℓ ≤ 20 range. Because of the latter, the largescale component of the field varies over the simulated sky, which is essential to reproduce the statistics of p. The other input sky components are taken from the latest version of the Planck Sky Model (Planck Collaboration XII 2016). These model component maps are then combined with the first 100 realizations of the systematic effects and noise (Planck Collaboration III 2020). This results in 100 E2E I, Q, and U maps at 353 GHz, which we smooth to 60′ resolution, and from which we derive polarization quantities and compare them to the input dust maps, after subtraction of the CMB and the CIB monopole, as is done for the ASMs (Sect. 2.4). This allows us to assess the effects of residual systematics and data noise on the statistics presented in this paper.
Figure C.1 shows that the difference between the input dust polarization fraction p^{(0)} and the average ⟨p⟩ over the 100 realizations of the E2E simulations is at most around 1%. The distribution function of the difference p − p^{(0)} is peaked around zero for each simulation. The average of these distributions is shown in the bottom panel of Fig. C.1, along with the 1σ dispersion around the average distribution, which has a mean of 0.03% and a standard deviation of 0.47%.
Fig. C.1. Top: map of polarization fractions p^{(0)} for the input sky of the E2E simulations, using the MAS estimator at 60′ resolution. Middle: map of the difference between the polarization fraction averaged over the 100 realizations of the E2E simulations, ⟨p⟩, and the input polarization fraction p^{(0)}, at 60′ resolution. Bottom: distribution function over the sky of the difference between the output and input polarization fractions. The solid blue curve is the average of 100 histograms of p − p^{(0)} from the 100 realizations, while the dashed lines with blue shading between show the ±1σ dispersion. 
For the polarization angle data, the diagnostic of the E2E simulations is shown in Fig. C.2. To compute the average difference between the output and input polarization angles, ⟨ψ⟩−ψ^{(0)}, account is taken of the circularity of the difference for each sky pixel and each realization independently. One can see that the regions of the sky where this difference is the largest are those where ψ^{(0)} crosses the ±90° boundary. The average distribution function of these differences over the 100 simulations and the 1σ dispersion about the average are shown in the bottom panel of the figure. The average distribution has a mean of 03 and a standard deviation of 83.
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. 
The same diagnostics are run on the polarization angle dispersion function 𝒮, which at 60′ resolution we compute with a lag δ = 30′. Results are shown in Fig. C.3. We note that the average ⟨𝒮⟩ of the output 𝒮 maps exhibits a significant positive bias with respect to the input 𝒮^{(0)} map, especially towards the Galactic poles. The distributions of pixel values in difference maps 𝒮 − 𝒮^{(0)} from the 100 simulations have a positive skewness. This shows that the polarization angle dispersion function is still affected by residual bias at this resolution, even though it is barely affecting the polarization angle map itself.
Fig. C.3. Top: map of the polarization angle dispersion function 𝒮^{(0)} for the input sky of the E2E simulations, at 60′ resolution and a lag of δ = 30′. Middle: map of the difference between the average ⟨𝒮⟩ and the input 𝒮^{(0)} in the E2E simulations. Bottom: distribution function of the difference between the output and input polarization angle dispersion function. The solid blue curve shows the average of 100 histograms of 𝒮 − 𝒮^{(0)} from the 100 simulations, with the 1σ dispersion shown as the blue area between dashed lines (barely visible). 
For completeness, Fig. C.4 shows the histograms of the polarization fractions, polarization angles, and polarization angle dispersions for the input (black curve), and the outputs of the E2E simulations. For the latter, the blue curve on each panel shows the average histogram over the 100 simulations, with the ±1σ dispersion among histograms shown as the blue area between dashed lines. The agreement is excellent for both quantities p and ψ, but for 𝒮 we note that at intermediate values the positive bias already mentioned appears clearly. Finally, we stress that although the polarization fractions rarely go above 20% for these simulated dust maps, this does not mean that the same range is expected in the Planck data.
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. 
Appendix D: Link between 𝒮 and the polarization gradients
D.1. Analytical derivation
The link between the polarization angle dispersion function 𝒮 and the polarization angle gradient ∇ψ can be established analytically via a Taylor expansion of the polarization angle difference appearing in the definition of 𝒮:
where the average is computed over the annulus centred on r having inner and outer radii δ/2 and 3δ/2, respectively (Eq. (4)), and ∇ψ is the vector gradient of the polarization angle at the centre r. Using a local reference frame with axes y and z in the plane of the sky, we can write the displacement vector as
with δ/2 ≤ l ≤ 3δ/2. The expression of 𝒮^{2} therefore becomes
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 cross 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
which leads to the relation
D.2. The case of Planck data
Figure D.1 shows the maps of both polarization gradients, ∇ψ and ∇P from Eqs. (6) and (5), respectively, for the GNILCprocessed Planck data at 353 GHz and 160′ resolution. The correlations between ∇ψ and 𝒮 on the one hand, and between ∇P and 𝒮 on the other, are shown in Fig. D.2 (for 𝒮 a lag of 80′ is used). These plots show that 𝒮 correlates well with ∇ψ, but not as well with ∇P and that ∇ψ is a very good proxy for the angular dispersion function 𝒮, as would be expected from Eq. (D.7) (and much faster to compute in practice).
Fig. D.1. Maps of the angular polarization gradient ∇ψ (left) and of the polarization gradient ∇P (right), built from the GNILCprocessed Planck data at 353 GHz and 160′ resolution. 
Fig. D.2. Left: twodimensional histogram representation of the correlation plot between the angular polarization gradient ∇ψ and the angular dispersion function 𝒮 from the GNILCprocessed Planck data at 353 GHz and 160′ resolution, with a lag of 80′ for 𝒮. Right: correlation plot between the polarization gradient ∇P and the angular dispersion function 𝒮. In both plots, the solid black curve shows the mean ∇ψ or ∇P in a given bin of 𝒮. 
Appendix E: Polarization fraction versus total gas column density for low and high offsets
In this appendix, we show in Fig. E.1 plots similar to Fig. 9, but for the low and high total intensity offsets (Sect. 2.2). The effects of the offset are discussed in Sect. 4.2.1.
Fig. E.1. Twodimensional 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. 
Appendix F: The inverse relationship between 𝒮 and p
In this appendix, we use a phenomenological model of the submillimetre polarized thermal dust emission, developed in Planck Collaboration Int. XLIV (2016), Ghosh et al. (2017), and Vansyngel et al. (2017), to derive the relationship between the polarization fraction p and the polarization angle dispersion function 𝒮. In its most basic form presented in Fig. F.1, this model assumes the polarized sky to result from a small set of N concentric layers, each emitting a fraction 1/N of the total intensity^{21}, and harbouring a magnetic field B = B_{0} + B_{t}, where B_{0} is a uniform field (the same in each layer) and B_{t} is an isotropic turbulent field that is taken, in each layer, as a different realization of a Gaussian random field in three dimensions. No effects of dust evolution or changes of intrinsic polarization properties of the dust grains are included in the model. By design, this model is able to reproduce the 1point statistics of polarized thermal Galactic dust emission maps observed by Planck, but it turns out that it is also able to reproduce the trend 𝒮 ∝ 1/p and the probability density function of 𝒮 × p, as we demonstrate below.
Fig. F.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. 
F.1. Reference frame and notations
We use a reference frame defined in Fig. F.2. The x axis is the line of sight, oriented towards the observer, and Oyz is the plane of the sky. In that frame, the components of the largescale magnetic field are (B_{0x}, B_{0y}, B_{0z}).
Fig. F.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. 
F.2. Magnetic field in a layer at a given line of sight
We begin by noting that, from one layer to the next, the different Gaussian realizations of the turbulent magnetic field can be taken to be independent. Therefore, in each layer i (with 1 ≤ i ≤ N), we can write the components of the magnetic field B_{i} = B_{0} + B_{t}, at the position considered as the central pixel in the definition of 𝒮 (Eq. (4)), as
Here G_{x}, G_{y}, and G_{z} are Gaussian random variables with zero mean and variance 1/3. The parameter
is the ratio of the standard deviation of the turbulent magnetic field to the magnitude B_{0} = B_{0} of the ordered field. The orientation of the magnetic field at the central pixel in layer i is given by a set of two angles Γ_{i} ∈ [ − π/2, π/2] and ϕ_{i} ∈ [0, 2π]:
As presented in Fig. F.2, the angle Γ_{i} is the inclination angle of the vector B_{i} with respect to the plane of the sky, while ϕ_{i} is the angle, counted positively clockwise from the north, between the z axis and the projection B_{POS} of B_{i} onto the plane of the sky.
F.3. Fluctuations within each layer over the scale δ
When computing the polarization angle dispersion function 𝒮, we introduce a specific scale, the lag δ, which we always take as half the FWHM ω, so that δ = ω/2. Presumably, the orientation of the magnetic field in each layer, i.e., the angles Γ_{i} and ϕ_{i}, vary little over these scales. Let us therefore consider a small Gaussian fluctuation δB_{i} around the direction of B_{i}. The rms σ_{Bi}(δ) of this fluctuation can be cast into a parameter similar in form to Eq. (F.4),
which depends on the lag δ considered, and is related to the overall turbulence parameter f_{M} and to the spectral index α_{M} of the magnetic field^{22}. This fluctuation δB_{i} corresponds to small variations for angles δΓ_{i} and δϕ_{i}:
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 is still small, which only fails if the direction of the mean magnetic field is very close to the line of sight.
F.4. Polarization angle and Stokes parameters
The angle ϕ_{i} is related to the polarization angle ψ_{i}, appearing in the definition of the Stokes parameters below, by a π/2 rotation, i.e., ψ_{i} = ϕ_{i} − π/2[π]. The πambiguity arises because the Stokes parameters are unchanged in the transformation B_{POS} ↦ −B_{POS}. The polarization angle thus lies in the range [ − π/2, π/2], and the Stokes parameters (I_{i}, Q_{i}, U_{i}) at the central pixel for each layer i are then^{23}
where I_{i} and P_{i} are the total and polarized intensity at the central pixel in layer i, respectively, and p_{max} is the polarization fraction of thermal dust emission that would be observed in the case of a uniform magnetic field parallel to the plane of the sky (Γ_{i} = 0). A fluctuation δB_{i} of the magnetic field at the scale δ therefore produces a small variation of these Stokes parameters that can be written as
where it is assumed that the total intensity varies little on the scale δ. Because we work with a lag smaller than the FWHM, δ = ω/2, this is a reasonable assumption. The fluctuation of the polarization angle is δψ_{i} = δϕ_{i}, and so by inserting Eqs. (F.12) and (F.13) for the fluctuations of the angles we obtain
These expressions will be helpful in determining the fluctuations of the Stokes parameters over which to average when computing the polarization angle dispersion function in the next section.
F.5. Polarization angle dispersion function
The polarization angle dispersion function 𝒮 is computed for a central pixel c, and consists of an average over the n pixels, indexed by j (with 1 ≤ j ≤ n), in an annulus of mean radius δ = δ and width δ around the central pixel, as defined in Eq. (4). This can also be written in terms of the Stokes parameters Q and U at the central pixel, and Q(j) and U(j) at a pixel j in the annulus (Planck Collaboration Int. XIX 2015):
Because we are interested in the average behaviour of 𝒮, we will ultimately consider the mean of this expression over the position of the central pixel as well.
The distribution function of 𝒮 (Fig. 7) shows that most pixels have a small dispersion of polarization angles, 𝒮 ≲ 10°. For these values, it is safe to approximate the arctangent by its argument, so that we may write
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
exhibiting the fluctuations of the Stokes parameters given in Eqs. (F.18) and (F.19). We use this decomposition to write the numerator and denominator that appear in the squared quantity above as
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
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}
Because the random variables g_{x}, g_{y}, and g_{z} are also uncorrelated from one layer to the next, we have
which yields, using the expressions of Eqs. (F.18) and (F.19),
Combining the above expressions, we then have
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:
We then have
which can be simplified further, using , to
Furthermore, in our phenomenological model the total intensity is split equally among the N layers, so that I_{i} = I/N. Therefore,
where 𝒜 is a geometrical factor that depends on the magneticfield structure in the layers, with
F.6. Application to Planck data: the case for strong turbulence
In line with our analysis of the data, we compute the mean of 𝒮(δ) over those pixels that have the same polarization fraction p. This gives
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 . In that case, the polarization angle dispersion function simply reads
We reach the important conclusion that the trend 𝒮 ∝ 1/p observed in the Planck data can be reproduced as a generic behaviour that depends only on the statistical properties of the turbulent magnetic field, without invoking changes in properties of the dust or in its alignment with respect to the magnetic field.
We note that the typical value and dispersion of the product 𝒮 × p depend not only on the properties of the turbulence at the scale of the lag, via the parameter, and on the number of layers N, but also on the maximum polarization fraction p_{max} that the dust can produce. Estimates of the latter are quite sensitive to the uncertainty on the zero level of the total intensity, as discussed in the main text.
For completeness, Fig. F.3 presents the probability density function (PDF) and cumulative density function (CDF) of various distribution functions. The solid curves correspond to the distribution of the meannormalized product 𝒮 × p/⟨𝒮×p⟩ for our model taken at a resolution of 160′ (models with up to 20′ resolution are similar). Empirically, the corresponding density functions for a Gamma distribution (Hazewinkel 2013) with shape parameter k = 5 and scale parameter θ = 1/5 reproduce these curves well, i.e., this Gamma distribution has similar statistics. These model and empirical density functions are also in reasonable agreement with the PDF and CDF of 𝒮 × p/⟨𝒮×p⟩ for Planck data at the same 160′ resolution.
Fig. F.3. Probability density function (PDF) and cumulative density function (CDF) of various distribution functions. Shown are the meannormalized product 𝒮 × p/⟨𝒮×p⟩ for our model taken at a resolution of 160′ (solid red curves), the same for Planck data at 160′ resolution (dot dashed orange curves), and a Gamma distribution with shape parameter k = 5 and scale parameter θ = 1/5 (dashed blue curves). 
F.7. Derivation of the expression for f_{m}(δ)
In our model, each component B_{ix}, B_{iy}, and B_{iz} of the magnetic field vector in layer i is the sum of the uniform field and a realization of a Gaussian random variable on the sphere, with a powerspectrum C_{ℓ} = C ℓ^{ αM}, where ℓ is the multipole. As in Planck Collaboration Int. XLIV (2016), we consider that the nonvanishing 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. (F.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 , 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:
where δB_{ix, ω} = B_{ix, ω} − B_{0x}, and similarly for the other components. The loss of power at large ℓ associated with the smoothing is clearly seen in Fig. F.4.
Fig. F.4. Power spectrum ℓ(ℓ + 1/2) 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 of the original power in the turbulent component (see text in Appendix F.8 and Eq. (F.54)). Shown as hatched is the turbulent energy implied in the calculation of 𝒮 at a resolution of ω = 80′ (with δ = ω/2 = 40′), which is a fraction of the original power in the turbulent component (see text in Appendix F.7 and Eq. (F.51)). Both coloured and hatched regions scale with the resolution ω as ω^{−2 − αM}. 
The factor f_{m}(δ) appearing in the expression of 𝒮 × p (Eq. (F.45)) is by definition (Eq. (F.8)) the typical relative fluctuation of the magnetic field at those scales comprised in the annulus between δ/2 and 3δ/2 (with δ = ω/2), i.e.,
Through Parseval’s theorem, and using the Fourier transform of this annulus of mean radius δ = δ and width δ (Gautier et al. 1992), we find (Fig. F.4)
where H(x) = 2J_{1}(x)/x, with J_{1} the Bessel function of the first kind of order one, and δ expressed in radians. Using Eq. (F.51), we compute at a resolution ω = 80′ (corresponding to δ = 40′) for α_{M} = −2.5. We find , which corresponds graphically to the hatched region in Fig. F.4.
To determine the dependence of on the lag δ, we simplify Eq. (F.51) by considering a unique ℓ ∼ 1/δ and a constant C_{ℓ} = C_{1/δ} in the sum, i.e., we replace the smooth, Gaussian, cutoff by a step in ℓ. This simplification, while being numerically incorrect, conserves the scaling of the integral with δ as long as ℓ ≫ 1.
Thus
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 .
Recalling δ = ω/2 to convert to ω and renormalizing to 160′, this analysis yields the following scaling with resolution:
valid as long as α_{M} does not depart too much from −2.5.
F.8. Beam depolarization
In this section, we estimate the effect of the resolution on the polarization fraction by quantifying the depolarization that occurs within the beam. This is important not only for comparing our results at 80′ and 160′ but also for taking into account the effects of the difference in resolution between the Planck polarization data and the starlight polarization that occurs within a pencil beam.
Following our approach in the previous section, we compute the difference in the total energy of the turbulent component, , between two given resolutions ω_{1} and ω_{2} > ω_{1}, for a given line of sight. We have
where σ_{1} and σ_{2} are related to ω_{1} and ω_{2} by , as already mentioned. We compute for ω_{1} = 80′, ω_{2} = 2ω_{1} = 160′ and α_{M} = −2.5. We find (see also Fig. F.4). Following the same approach as for Eq. (F.52), this yields the following scaling with the resolution ω:
From Eqs. (F.53) and (F.55), we conclude that the factor is independent of f_{M} and ω, and only depends on α_{M}. For α_{M} = −2.5, it is equal to 3.02.
We now study the effect of smoothing on the polarization fraction map. For that we note that Stokes Q and U exhibit a power spectrum similar to that of the turbulent component of the magnetic field at ℓ ≫ 1 (see Appendix F.9). From Parseval’s theorem, we therefore have, for each layer i
where the factor comes from the loss of power in the turbulent component of the field between the two resolutions (Fig. F.4) and k is a constant to be determined numerically. The different random realizations δB_{i} are independent from one another for ℓ ≥ 2. At small enough spatial scales (ℓ ≫ 1), this ensures that the various realizations of Q_{i} and U_{i} are also independent from one another. We therefore have
A comparison with our numerical model shows that k ≃ 1 (Fig. F.5 presents the agreement with the model for p when taking k = 1), and therefore
Fig. F.5. 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 in the rms of p (red) by depolarization (Eq. (F.62) with ω = 160′) and the increase in 𝒮 × p (black) with the resolution (Eqs. (F.53) and (F.45)). The fractional difference is less than 10% for 𝒮 × p. The dashed blue line represents the pencil beam value for the rms of p, as calculated from Eq. (F.63) based on the model taken at ω = 160′. 
Using Eq. (F.45), this yields
Equation (F.59) quantifies by how much the polarization fraction decreases when smoothing from resolution ω to 2ω, because on average (only) .
We can now generalize to the case of smoothing data from a finer resolution ω/2^{n} to a resolution ω. In that case, we can compute the beam depolarization by a chain sum,
where all averages are taken over the entire map. This yields
In summary, the change in squared polarization fraction from the coarser scale ω to the finer resolution ω/2^{n} is
In the pencilbeam limit (n = ∞), we have
From symmetry arguments, we also have
and the same for U/I.
F.9. Comparison of the analytical expressions with numerical results and application to pencil beams
In Fig. F.5, we compare our analytical expressions for the mean 𝒮 × p (Eqs. (F.53) and (F.45)) and the rms of p (Eq. (F.62)) as a function of the resolution, with numerical results directly computed from the smoothed Stokes maps of our simulated model from Planck Collaboration Int. XLIV (2016), i.e., f_{M} = 0.9, α_{M} = −2.5, N = 7 and p_{0} = 26%. There are two aspects of the comparison for the analytical model, the normalization and the dependence on resolution.
For 𝒮 × p, we observe a slight normalization difference (7%) between the analytical and numerical results. More precisely, at 160′, we find 𝒮 × p = 034 for the analytical expression and 032 for the simulation. Both are nevertheless very close to the observational value of 031. A deviation from the trend predicted by Eq. (F.53) is observed at low resolutions, as expected from our demonstration in Appendix F.7.
Concerning the beam depolarization, the decrease in the polarization fraction with a decreasing resolution (larger ω) is well reproduced over two orders of magnitude in resolution: the example shown corresponds to a mean rms of p of 11% over the full sky. Nevertheless, Vansyngel et al. (2017) already noted a small (approximately 0.1) difference between the index α_{M} characterizing the power spectrum of the turbulent component of the magnetic field in the simulation, and the index α_{EE} and α_{BB} recovered from the analysis of the EE and BB power spectra. This is also what we find here: the scaling of f_{m} and f_{p} with ω is actually closer to ω^{0.18}, which would correspond to α_{M} = −2.36, when the model is produced with α_{M} = −2.5. We show this scaling as the solid lines in Fig. F.5 and this is why in the rest of the paper we consider a scaling
The case of a pencil beam is of interest. Applying Eq. (F.63) to the highest polarization fraction observed at 160′, p_{max}(ω = 160′)≈20%, we can estimate the corresponding rms polarization over that scale in pencilbeam data, which corresponds to n = ∞ and for which the prefactor in Eq. (F.63) is 63.9 for α_{M} = −2.36. Expressing the observed ⟨𝒮 × p⟩_{160′} = 031 in radians, we find an rms p_{max}(ω = 0′) = 20.5%. In that particular case, depolarization is expected to be very small because when p = p_{max}, the magnetic field is already ordered and within the plane of the sky and therefore little subject to depolarization. However, the effect would be stronger for lines of sight characterized by a low polarization fraction. For example, if p(ω = 160′) = 6.0% at 160′ resolution, Eq. (F.63) gives an rms polarization fraction p(ω = 0′) = 7.4% over that scale for pencilbeam data.
Appendix G: Noise and systematics in the comparison of submillimetre and optical polarization data
For our estimation of the emissiontoextinction polarization ratios R_{P/p} and R_{S/V} (Sect. 6), all observations should ideally be done in a pencil beam and probe the full line of sight through the Galactic dust. Unavoidable departures from this ideal situation therefore introduce systematic effects on the quantities appearing in these ratios. In this appendix, we estimate these various systematic effects.
G.1. Beam depolarization
Systematic distortions of the submillimetre polarization signal occur due to the averaging of Stokes Q and U in the telescope beam. This does not happen with the pencil beam of optical observations. In Appendix F.8 we demonstrate that beam depolarization produces a negative bias in p with respect to the pencilbeam value, with a scatter around this biased value (Eq. (F.63)). For α_{M} = −2.36, we find (Appendix F.9). As a consequence, to compare optical and submillimetre polarization data, we compensate for this systematic beam depolarization by multiplying all Planck Stokes parameters and uncertainties (P, p, σ_{P}, σ_{p}) by the factor . Because Q and U play a symmetric role in P^{2} = Q^{2} + U^{2}, the same correction factor is applied to Q, U, Q/I, and U/I.
G.2. Background distortion
As mentioned in Sect. 6.2, the optical polarization degree p_{V} is potentially biased by the difference in length probed along the line of sight, compared to the polarized emission in the submillimetre. Under the assumption of a uniform ISM, this bias due to the background may be corrected using Eq. (18). Recognizing that the background is not uniform, we estimate that the uncertainty on (the ideal measure from the infinity) is proportional to the amount of reddening behind the star ΔE(B − V)^{∞, ⋆} = E(B − V)^{∞} − E(B − V)^{⋆} and to the rms of the polarization degree per unit reddening in the background, i.e.,
This equation replaces Eq. (19) when the background is not uniform. We discuss the appropriate value of the last factor in Appendix G.4.1.
For the polarization fraction p_{V}/τ_{V}, no renormalization is needed as both p_{V} and τ_{V} are measured from the star. The systematic uncertainty on p_{V}/τ_{V} from background distortion in then simply taken as:
These equations also apply respectively to uncertainties on and (Eq. (G.1)), and q_{V}/τ_{V} and u_{V}/τ_{V} (Eq. (G.2)).
G.3. Uncertainties related to the reddening maps
The uncertainty on E(B − V)^{⋆} stems from the uncertainty on the PS1based reddening data and the uncertainty σ_{θ⋆} on the stellar parallax. We estimate the former by correlating the PS1 total reddening, E(B − V)^{∞}, with the Planck optical depth at 353 GHz converted to a reddening, E(B − V)_{τ}. Figure G.1 shows that these quantities are remarkably well correlated, with a Pearson correlation coefficient r = 0.98 and little scatter (around 15% of the running mean on average). Lacking more precise information, we assume equal contributions from instrinsic scatter between E(B − V)^{∞} and E(B − V)_{τ} on the one hand and from noise in E(B − V)^{∞} on the other hand, so that we take σ_{E(B − V)}/E(B − V) = 0.1 for the PS1based estimates. The small departure from linearity in this correlation has been interpreted as evidence for dust evolution in the diffuse ISM (Planck Collaboration XI 2014; Planck Collaboration Int. XXIX 2016), i.e., for the modification of the dust optical properties in its lifecycle through the ISM. This aspect will be investigated in its relation to dust polarization properties in a future paper.
Fig. G.1. 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, for the 1505 selected stars (Sect. 6.3). Each bin of the running mean contains the same number of lines of sight. Error bars represent the standard deviation in each bin. The dashed line corresponds to a onetoone correlation. 
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 and are the reddenings to the star obtained for the altered parallaxes θ^{⋆} − σ_{θ⋆} and θ^{⋆} + σ_{θ⋆}, respectively. Gathering the two sources of uncertainty, the total uncertainty on E(B − V)^{⋆} is then
This uncertainty then propagates to the quantities used in determining R_{S/V} (see Sect. 6.4), e.g.,
and similarly for q_{V}/τ_{V} and u_{V}/τ_{V}.
G.4. Polarization angle difference
The unbiased comparison between submillimetre and optical measurements also requires an agreement in polarization angles. We recall that for interstellar polarization of starlight, the direction of the projection of the magnetic field on the plane of the sky (B_{POS}) can be inferred directly from the polarization angle. For polarization in emission at 353 GHz, a rotation by 90° is required. In Fig. G.2, we compare the direction of B_{POS} inferred from the two tracers. The length of each line is proportional to the S/N value for the corresponding polarization angle. If the same dust is probed in the optical and the submillimetre, the directions should be identical, and we do find that the agreement is quite remarkable, although not perfect.
Fig. G.2. 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. 
To quantify this agreement, we define the difference in orientation angles between the Planck (or submillimetre, “S") and optical (or visual, “V") polarization data as Δψ_{S/V} = (ψ+90°) − ψ_{V}. In terms of Stokes parameters this can be written as (Planck Collaboration Int. XXI 2015)
Figure G.3 presents the histogram of Δψ_{S/V}, of which two important aspects should be understood, i.e., its width and any shift relative to the expected centring on zero.
Fig. G.3. Histogram of the difference in polarization angles between Planckderived angles and opticalpolarizationderived angles, Δψ_{S/V} = (ψ+90°) − ψ_{V}, with its standard deviation σ_{S/V} and median value indicated. Histograms are overplotted for simulations based on noise only (dashed red curve), and noise plus systematics (dashed blue curve) – see text. For easier comparison, the latter histogram has been shifted horizontally by −31 to peak at the same position as the data. 
G.4.1. Standard deviation
We attempt to explain the width of the histogram Δψ_{S/V} (Fig. G.3) by simulating what is expected from noise and systematic effects. We start from Planck data at a resolution of 40′ and starlight data, rotated to assume perfect orthogonality. Then we add fluctuations drawn from Gaussian distributions having the estimated variances for each random and systematic uncertainty considered among the following: noise in the submillimetre and optical, beam depolarization (Appendix G.1), and background distortion (Appendix G.2). Finally, with the resulting Stokes parameters, Q and U at 353 GHz, and q_{V} and uv in the optical, we compute Δψ_{S/V} for these simulated data.
The width of the observed histogram is too large to be explained by noise only, as shown by the dashed red curve that results from simulating only effects of the noise in the submillimetre and optical data (Fig. G.3). Part of this discrepancy may be accounted for by the dispersion of polarization angles in the optical within a Planck beam, i.e., at scales that cannot be probed by Planck. To estimate this contribution, we analyse the starlight polarization data alone, ignoring Planck data. We compute the standard deviation σ_{V/V} of the difference in polarization angles in the optical, for those stars at an angular distance δ smaller than 40′ (one FWHM of the Planck beam). We find σ_{V/V} = 26°. We note that by its nature this dispersion incorporates twice the variance of optical polarization angles compared to that which enters into the Δψ_{S/V} histogram. Thus the expected standard deviation from this effect would be about 18°. Three random and systematic effects affecting starlight polarization measures could explain that startostar dispersion : noise in the optical, turbulence at scales smaller than the Planck beam (Appendix G.1); and background distortion (Appendix G.2). For these three, we consider Gaussian fluctuations having variances , ^{24} and , respectively, the latter depending on an unknown parameter, rms(p_{V}/E(B − V))_{bkgd}. We find that simulations including both optical noise and turbulence within the Planck beam produce a startostar dispersion of 13°, well below the observed value of 18°. This points to the need for a contribution from background distortion, and thus a need to estimate rms(p_{V}/E(B − V))_{bkgd}. The present data sample provides lower and upper bounds to the rms of p_{V}/E(B − V), respectively 7% for the rms of p_{V}/E(B − V) in our sample of 1505 stars, and 13% for the maximum observed polarization fraction (see Sect. 6.6). Assuming a value of rms(p_{V}/E(B − V))_{bkgd} = 8% close to the rms over our sample of stars, we obtain a dispersion of 126 arising from the background distortion . The total dispersion including optical noise, turbulence within the Planck beam, and background distortion is then 16°, still significantly below our estimate of 18°. We note that the individual contributions to the angular dispersion do not quite add up in quadrature to the one obtained when taking into account all sources of errors, because even though the fluctuations of Stokes parameters are Gaussiandistributed, the polarization angles are not (NaghizadehKhouei & Clarke 1993). A total dispersion of 18° would require rms(p_{V}/E(B − V)_{bkgd}) = 11%, a choice that is probably too extreme given the observational constraints, the uncertainties on our estimated value for σ_{V/V}, and other physical effects not included in this analysis (see Appendix G.4.2 below). Therefore, to be conservative, we will adopt a value of 8% for the rms of p_{V}/E(B − V)_{bkgd} to compute the background distortion, which then by itself contributes a dispersion of 126.
Our simulated histogram of Δψ_{S/V} based on optical and submillimetre noise, beam depolarization at 40′ and background distortion is presented in Fig. G.3 (dashed blue line). It is close in shape and width (standard deviation σ_{N + S} = 208) to the observed histogram and by construction it is centred near zero. The contributions to the standard deviation are 139 from submillimetre noise, 88 from optical noise, 96 from beam depolarization, and 126 from background distortion. Here again, we warn against the simple quadratic addition of individual sources of uncertainty. We note that this model does not include a contribution from possible Planck systematics, which can be assessed through the E2E simulations presented in Appendix C. For the lines of sight to stars, we built the histogram of the difference in polarization angles between the input maps at 5′ resolution and those at 40′ that went through the simulation pipeline and therefore include estimates of Planck noise and known systematics. The standard deviation of these histograms is 14° ± 1°, depending on the number of simulations used, a value close to that found from Planck noise alone, underlining the low level of remaining systematics in the Planck data.
This analysis of the standard deviation of Δψ_{S/V} provides some assurance that the subset of uncertainties needed for the quantification of the emissiontoextinction polarization ratios in Sect. 6.4 have been quantified adequately.
G.4.2. Mean difference
The histogram of Δψ_{S/V} in the data peaks at −31, revealing a systematic angle difference (shift) between polarization angle measurements in emission relative to extinction. Given the large number of lines of sight, this shift cannot be explained by a random process, thus pointing to a systematic effect. Although the shift is small and unimportant for evaluating the polarization ratios discussed in Sect. 6, it is potentially important for other investigations. For example, accurate absolute calibration of the polarization angle is critical for future CMB Bmode experiments to avoid systematic effects that could compromise reaching the precision required (Aumont et al. 2018). We have explored the possible origins of such a shift.
There is some evidence for the possibility of background distortion, a systematic difference arising because in the submillimetre dust is observed along a longer path length than probed in the optical. First, within the full 2388star sample the shift appears to depend on Galactic longitude. Second, when we apply a more stringent criterion on the reddening ratio E(B − V)^{⋆}/E(B − V)^{∞} (see Sect. 6.3), thus reducing the chance of a significant background contamination, we find that the shift is smaller, −2°, and a longitude dependence is no longer evident.
In the abovementioned Planck E2E simulations (Appendix G.4.1), the mean and median of the difference in polarization angles are consistently shifted from zero by less than 025, which shows that the observed difference is not ascribable to any known systematic in the Planck data. Concerning the calibration of the zero point of the polarization angle in the submillimetre, the uncertainties on the orientation of the HFI PSBs at 100, 143, and 217 GHz are below 03 (Planck Collaboration Int. XLVI 2016; Planck Collaboration XI 2020), and it seems likely that the uncertainty at 353 GHz is of the same magnitude. This is consistent with multifrequency measurements of the polarization angle of the bright synchrotron emission of the Crab Nebula from 23 to 353 GHz (Aumont et al. 2018).
The calibration of the zero point in the optical is less clear. For the highlatitude polarization surveys considered here, three highly polarized stars were used (Berdyugin et al. 2014). We have examined the extensive historical record of the polarization angle of these stars and find good evidence that each varies with an rms of typically 15, but with excursions as large as 5°. Without confidence in subdegree accuracy of the optical zero point, the possibility remains of a significant contribution to the shift due to this uncertain calibration. There is not a lot of overlap between the optical polarization compilation of Heiles (2000) and the surveys of Berdyugin et al. (2014) used here. However, by a statistical comparison between the polarization angles of stars in one catalogue with those of stars in the other, binned as a function of angular distance, it is possible to investigate any systematic shift. This analysis can be extended by replacing the optical polarization angle measurements of one or both catalogues with the Planck polarization angle at the catalogue positions. This reveals small systematic shifts up to a few degrees that could arise from different zero point calibrations and/or different path lengths probed. From these investigations, a residual contribution to the shift of order 1° seems possible.
A final possibility is actual decorrelation between the submillimetre and optical polarization, e.g., due to a temperatureweighting effect in the submillimetre coupled with a correlation of variations in the heating of aligned grains (or in the dust properties themselves) along the line of sight with variations of the magnetic field orientation (Tassis & Pavlidou 2015).
We cannot pinpoint a single cause for a shift of the magnitude seen. Instead, it seems that several smaller contributions might have conspired to produce the effect observed.
All Tables
All Figures
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. The scales for I and P are logarithmic, while those for Q and U are linear. 

In the text 
Fig. 2. Polarization maps for the GNILCprocessed 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. 

In the text 
Fig. 3. Signaltonoise ratio (S/N) p/σ_{p} for the polarization fraction in the GNILCprocessed data at 353 GHz and uniform 80′ resolution. The polar view (bottom) uses a range 1 ≤ p/σ_{p} ≤ 10 to bring out low S/N regions. 

In the text 
Fig. 4. Top: polarization angle dispersion function 𝒮 computed from the GNILCprocessed data at 353 GHz and uniform 160′ resolution, using a lag δ = 80′. Bottom: statistical uncertainty computed from the Monte Carlo simulations on maps with the same 160′ resolution and δ = 80′ lag. 

In the text 
Fig. 5. Distribution functions of the polarization fraction p in the GNILC data at 353 GHz and uniform 80′ resolution. The solid red 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. 

In the text 
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 data, the dashed curve gives the mean of the 1000 Monte Carlo histograms, and the blue region shows the envelope spanned by the 1000 histograms. 

In the text 
Fig. 7. Distribution functions of the polarization angle dispersion function 𝒮 in the GNILC data at 353 GHz. The cases shown are for the 160′ resolution using a lag δ = 80′ (in green), and for the 80′ resolution using a lag δ = 40′ (in blue). The solid curves show the histograms computed directly from the GNILC maps, the dashed curves give the mean histogram from the 1000 Monte Carlo realizations for each case, and the coloured regions show the envelope. The dashed vertical line indicates the value corresponding to a completely random polarization pattern. 

In the text 
Fig. 8. Distribution functions of 𝒮 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. 

In the text 
Fig. 9. Twodimensional 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. 

In the text 
Fig. 10. Twodimensional histogram showing the joint distribution function of 𝒮 and p at 160′ resolution, using a lag δ = 80′. The black curve is the running mean of 𝒮 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 𝒮 in each bin of p. The dashed white line shows our fit 𝒮 = 031/p to this running mean. 

In the text 
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. 

In the text 
Fig. 12. Polarization fraction p as a function of the total gas column density, N_{H}, coloured by the polarization angle dispersion function 𝒮 (on a logarithmic scale). The resolution of the data is 160′, in order to limit the bias in 𝒮. 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. 

In the text 
Fig. 13. Polarization fraction p (top panel), and polarization angle dispersion function 𝒮 (bottom panel) coloured by the dust temperature, T_{d}, as a function of the column density N_{H}. We note that the resolution of the data used here is 40′. Estimates of 𝒮 are nevertheless not biased, because only lines of sight with N_{H} ≥ 10^{21} cm^{−2} are considered. 

In the text 
Fig. 14. Twodimensional histograms with background colours encoding the density of points on a logarithmic scale, showing p (top), 𝒮 (middle), and 𝒮 × p (bottom) as a function of the column density N_{H} (left), Galactic latitude b (middle), and Galactic longitude l (right). The resolution is 160′. The colour bar shown in the top left panel is common to all plots. Black curves show the running means calculated as in Fig. 12, with error bars representing the scatter in each bin. For 𝒮, which is on a logarithmic scale, the median trend shown (thin black line) follows the density of points more faithfully than does the mean (thicker black line). 

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

In the text 
Fig. 16. Mean 𝒮 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 𝒮 = f(p) trend at 160′ (Fig. 10), corrected for its dependence on the resolution, as per Eq. (10). All bins in p contain the same number of pixels, n ≈ 250. Error bars correspond to the uncertainty on the mean, i.e., , where σ is the statistical dispersion in the corresponding bin. 

In the text 
Fig. 17. Correlation between dust temperature T_{d} and our estimate G_{ℛ} 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). 

In the text 
Fig. 18. Polarization fraction at 353 GHz (top) and product 𝒮 × 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 𝒮 × p at 160′ (Fig. 10), corrected for its dependence on the lag (Eq. (10)). 

In the text 
Fig. 19. Mean of p (top), 𝒮 (middle), and 𝒮 × 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 colour represents the density of points at a resolution of 10′, shown only for N_{H} > 4 × 10^{21} cm^{−2}. 

In the text 
Fig. 20. Mean S/N of p (top), and 𝒮 (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, not to the uncertainty on the mean. The dashed line indicates the minimal S/N that ensures reliable mean values for debiased quantities. 

In the text 
Fig. 21. Mean 𝒮, p, and 𝒮 × p as a function of N_{H}, combining results from Planck maps at optimal resolutions for all lines of sight above b > 5° (solid curves). For clarity, 𝒮 has been raised vertically by a factor of 2. Upper and lower dashed red 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 on the mean, , 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 thin black dashed (b > 10°) and dotted (b > 2°) curves. Horizontal coloured bars indicate for each resolution ω the column density interval ℐ(ω) used in the renormalization procedure (see text). The green band highlights a 25% decrease in 𝒮 × p with column density up to 2 × 10^{22} cm^{−2}. 

In the text 
Fig. 22. Histogram of the ratio of the reddening to the star to the total reddening on the same line of sight, E(B − V)^{⋆}/E(B − V)^{∞}, as derived from the PanSTARRS1 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. 

In the text 
Fig. 23. Correlation between Stokes polarization parameters in emission at 353 GHz and in optical extinction, with the colour in the 2dimensional histogram representing the density of points. Left: Stokes parameters (Q, U) versus (), 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). The reduced χ^{2}, the Pearson correlation coefficient, and the correlation coefficient inferred from the Bayesian method (Kelly 2007) are listed. 

In the text 
Fig. 24. Emissiontoextinction 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}. For the running mean each bin contains the same number (100) of lines of sight. The lower dark blue dotteddashed lines indicate the reduced χ^{2} of the fits (with the scale on the left axis). Dashed red and purple curves represent the Pearson and Bayesian (Kelly 2007) 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). 

In the text 
Fig. 25. Polarization fraction in emission at 353 GHz, p = P/I (left), and in optical extinction, p_{V}/E(B − V) (right), as a function of the column density, N_{H}. The sample in blue shows 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 (1505 stars), where data have been corrected for systematic effects such as beam depolarization (see Appendix G). 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) ∝ E(B − V)^{−0.2}, is shown for comparison. 

In the text 
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′. The scale of the I map is logarithmic, while the rest are linear. 

In the text 
Fig. B.2. GNILCprocessed covariance maps at 353 GHz and varying resolution. From top to bottom and left to right, they are σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU}. 

In the text 
Fig. B.3. GNILCprocessed covariance maps at 353 GHz and unifom 80′ resolution. From top to bottom and left to right, they are σ_{II}, σ_{IQ}, σ_{IU}, σ_{QQ}, σ_{QU}, and σ_{UU}. 

In the text 
Fig. C.1. Top: map of polarization fractions p^{(0)} for the input sky of the E2E simulations, using the MAS estimator at 60′ resolution. Middle: map of the difference between the polarization fraction averaged over the 100 realizations of the E2E simulations, ⟨p⟩, and the input polarization fraction p^{(0)}, at 60′ resolution. Bottom: distribution function over the sky of the difference between the output and input polarization fractions. The solid blue curve is the average of 100 histograms of p − p^{(0)} from the 100 realizations, while the dashed lines with blue shading between show the ±1σ dispersion. 

In the text 
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. 

In the text 
Fig. C.3. Top: map of the polarization angle dispersion function 𝒮^{(0)} for the input sky of the E2E simulations, at 60′ resolution and a lag of δ = 30′. Middle: map of the difference between the average ⟨𝒮⟩ and the input 𝒮^{(0)} in the E2E simulations. Bottom: distribution function of the difference between the output and input polarization angle dispersion function. The solid blue curve shows the average of 100 histograms of 𝒮 − 𝒮^{(0)} from the 100 simulations, with the 1σ dispersion shown as the blue area between dashed lines (barely visible). 

In the text 
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. 

In the text 
Fig. D.1. Maps of the angular polarization gradient ∇ψ (left) and of the polarization gradient ∇P (right), built from the GNILCprocessed Planck data at 353 GHz and 160′ resolution. 

In the text 
Fig. D.2. Left: twodimensional histogram representation of the correlation plot between the angular polarization gradient ∇ψ and the angular dispersion function 𝒮 from the GNILCprocessed Planck data at 353 GHz and 160′ resolution, with a lag of 80′ for 𝒮. Right: correlation plot between the polarization gradient ∇P and the angular dispersion function 𝒮. In both plots, the solid black curve shows the mean ∇ψ or ∇P in a given bin of 𝒮. 

In the text 
Fig. E.1. Twodimensional 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. 

In the text 
Fig. F.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. 

In the text 
Fig. F.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. 

In the text 
Fig. F.3. Probability density function (PDF) and cumulative density function (CDF) of various distribution functions. Shown are the meannormalized product 𝒮 × p/⟨𝒮×p⟩ for our model taken at a resolution of 160′ (solid red curves), the same for Planck data at 160′ resolution (dot dashed orange curves), and a Gamma distribution with shape parameter k = 5 and scale parameter θ = 1/5 (dashed blue curves). 

In the text 
Fig. F.4. Power spectrum ℓ(ℓ + 1/2) 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 of the original power in the turbulent component (see text in Appendix F.8 and Eq. (F.54)). Shown as hatched is the turbulent energy implied in the calculation of 𝒮 at a resolution of ω = 80′ (with δ = ω/2 = 40′), which is a fraction of the original power in the turbulent component (see text in Appendix F.7 and Eq. (F.51)). Both coloured and hatched regions scale with the resolution ω as ω^{−2 − αM}. 

In the text 
Fig. F.5. 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 in the rms of p (red) by depolarization (Eq. (F.62) with ω = 160′) and the increase in 𝒮 × p (black) with the resolution (Eqs. (F.53) and (F.45)). The fractional difference is less than 10% for 𝒮 × p. The dashed blue line represents the pencil beam value for the rms of p, as calculated from Eq. (F.63) based on the model taken at ω = 160′. 

In the text 
Fig. G.1. 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, for the 1505 selected stars (Sect. 6.3). Each bin of the running mean contains the same number of lines of sight. Error bars represent the standard deviation in each bin. The dashed line corresponds to a onetoone correlation. 

In the text 
Fig. G.2. 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. 

In the text 
Fig. G.3. Histogram of the difference in polarization angles between Planckderived angles and opticalpolarizationderived angles, Δψ_{S/V} = (ψ+90°) − ψ_{V}, with its standard deviation σ_{S/V} and median value indicated. Histograms are overplotted for simulations based on noise only (dashed red curve), and noise plus systematics (dashed blue curve) – see text. For easier comparison, the latter histogram has been shifted horizontally by −31 to peak at the same position as the data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.