© Planck Collaboration 2020

1. Introduction

The polarization of the cosmic microwave background (CMB) offers an opportunity for detecting primordial gravitational waves, a key experimental manifestation of quantum gravity (Starobinskiǐ 1979). Inflation generates tensor (gravitational waves) together with scalar (energy density) inhomogeneities. The polarization curl-like signal, referred to as primordial B modes, is a generic signature of gravitational waves produced during the inflation era in the very early Universe (Guth 1981; Linde 1982). However, the ratio of tensor-to-scalar power, denoted r, varies considerably among models (Baumann 2009). Improvement of the present limit, r <  0.07 (95% confidence, BICEP2 and Keck Array Collaborations 2016), might be achieved by combining data from new sub-orbital experiments with data from Planck¹, as pioneered by the BICEP/Keck and Planck joint analysis (BICEP2, Keck Array and Planck Collaborations 2015).

Until the next CMB space mission, the Planck data will remain unique, both for the all-sky coverage, required to measure CMB polarization at very low multipoles, and for its sensitive 353 GHz dust polarization maps. At microwave frequencies, the sensitivity of Planck is limited by the small number of detectors (12 per channel for the High Frequency Instrument HFI), while today the most sensitive sub-orbital experiments have array sizes up to of order 10³ detectors. Further in the future, the CMB stage III and IV development plans in the United States include array sizes increasing to more than 10⁵ detectors, with a goal of detecting primordial B modes down to r  ≃  10⁻³. On-going sub-orbital projects, including Advanced ACTPol (Naess et al. 2014), BICEP2/3 and the Keck Array (Grayson et al. 2016), CLASS (Essinger-Hileman et al. 2014), PIPER (Kogut et al. 2011), POLARBEAR and the Simons Array (Arnold et al. 2014), the Simons Observatory², SPIDER (Fraisse et al. 2013), and SPTPol (Austermann et al. 2012), are paving this ambitious path.

Indeed, the primordial B modes might have high enough amplitude to be discovered by these experiments, but this exciting prospect does not depend solely on the data sensitivity. Discovery depends on component separation, because the cosmological signal is contaminated by polarized foreground emission from the Galaxy that has a higher amplitude (Dunkley et al. 2009; BICEP2, Keck Array and Planck Collaborations 2015; Errard et al. 2016; Hensley & Bull 2018; Remazeilles et al. 2018). Component separation is also a key issue in the definition of future CMB space experiments, for example LiteBIRD (Ishino et al. 2016). This component-separation challenge binds the search for primordial B modes to the statistical characterization, and the astrophysics, of polarized emission from the magnetized interstellar medium (ISM).

The spin axis of a non-spherical dust grain is both perpendicular to its long axis and aligned, statistically, with the orientation of the ambient Galactic magnetic field. This alignment makes dust emission polarized perpendicular to the magnetic field projection on the plane of the sky (Stein 1966; Hildebrand 1988; Martin et al. 2007), and also perpendicular to the optical interstellar polarization from the same grains, as confirmed by Planck Collaboration Int. XXI (2015). Dust emission is the dominant polarized foreground at frequencies larger than around 70 GHz (Dunkley et al. 2009; Planck Collaboration X 2016). The Planck maps greatly supersede, in sensitivity and statistical power, the data available from earlier ground-based and balloon-borne observations.

Several studies have already used the Planck data to investigate the link between the dust polarization maps and the structure of the ISM and of the Galactic magnetic field (GMF). Planck Collaboration Int. XIX (2015) presented the first analysis of the polarized sky as seen at 353 GHz (the most sensitive Planck channel for polarized thermal dust emission), focusing on the statistics of the polarization fraction and angle, p and ψ. Comparison with synthetic polarized emission maps, computed from simulations of magneto-hydrodynamical (MHD) turbulence, shows that the turbulent structure of the GMF is able to reproduce the main statistical properties in interstellar clouds (Planck Collaboration Int. XX 2015).

Planck Collaboration Int. XXX (2016; hereafter PXXX) present the polarized dust angular power spectra computed with the Planck data over the high-Galactic-latitude sky that is best suited for the analysis of CMB anisotropies. An E/B asymmetry (usually quantified as the power ratio ) was discovered, as well as significant TE power. A correlation between the filamentary structure of cold gas identified in the Planck dust total intensity maps, and the local orientation of the GMF, derived from the dust polarization angle, has shown the two fields to be aligned statistically (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXV 2016; Kalberla et al. 2016). This alignment has also been reported for filamentary structures identified in spectroscopic HI data cubes (McClure-Griffiths et al. 2006; Clark et al. 2014, 2015). The structures identified in HI channel maps could, at least partly, correspond to gas velocity caustics. In that case, the correlation between gas velocity and magnetic field orientation (Lazarian et al. 2018) would contribute to the observed alignment. However, the Planck dust total intensity maps trace the dust column density, and for these data the observed correlation with the GMF is unambiguously an alignment of density structures with the magnetic field. Planck Collaboration Int. XXXVIII (2016) showed that this correlation could account for the E/B asymmetry and also the TE correlation.

These observational results have been discussed in the context of interstellar turbulence. The alignment between density structures and magnetic field is observed in MHD simulations of the diffuse ISM and discussed by Hennebelle (2013), Inoue & Inutsuka (2016) and Soler & Hennebelle (2017). The E/B asymmetry and the TE correlation have been considered as statistical signatures of turbulence in the magnetized ISM from different theoretical perspectives by Caldwell et al. (2017), Kandel et al. (2017, 2018), Kritsuk et al. (2018). This hypothesis is still debated. There is no consensus on whether it holds, and what we may be learning about interstellar turbulence.

Planck Collaboration Int. XLIV (2016) introduced a phenomenological framework that relates the dust polarization to the GMF structure, its mean orientation and a statistical description of its random (turbulent) component. This framework has been used to model dust polarization power spectra and to produce simulated maps that can be used to assess component-separation methods and residuals in the analysis of CMB polarization (Ghosh et al. 2017; Vansyngel et al. 2017) and also underlies the dust sky model in the end-to-end (E2E) simulations used in this paper (see Appendix A).

The Planck data on polarized thermal dust emission allowed Planck Collaboration Int. XXII (2015) to determine the spectral energy distributions (SEDs) of dust polarized emission and dust total intensity at microwave frequencies (ν ≤ 353 GHz). The combination of BLASTPol submillimetre data with Planck (Gandilo et al. 2016; Ashton et al. 2018) also shows that the frequency dependence of the polarization fraction p is not strong. New constraints like this, along with the ratio of dust polarized emission to the polarization fraction of optical interstellar polarization (Planck Collaboration Int. XXI 2015), can be used to refine dust models (Guillet et al. 2018). The modelling of the dust SED is also essential to component separation for CMB studies (Chluba et al. 2017; Hensley & Bull 2018).

Planck Collaboration Int. L (2017; hereafter PL) studied the correlation between dust polarization maps from the HFI channels at 217 and 353 GHz. In developing the analysis for this paper, we found that systematic errors and noise question the evidence for spectral decorrelation proposed in that earlier paper. This conclusion is in agreement with the results of Sheehy & Slosar (2018), who discovered this independently.

In this paper, one of a set associated with the 2018 release of data from the Planck mission (Planck Collaboration I 2020), we make use of Planck maps from this third public release hereafter referred to as PR3, to extend the characterization of the polarized Galactic dust emission that is foreground to CMB polarization. Our data analysis procedure has three main new directions.

• (1)

  We expand the power-spectrum analysis of dust polarization into the low-multipole regime relevant for E- and B-mode CMB polarization associated with the reionization of the Universe. This part of our analysis includes a validation of the dust polarization maps through running the mapmaking pipeline on simulated time-line data built from simulations of the sky, including a model of polarized dust emission.

• (2)

  We characterize the mean SED of polarized Galactic foregrounds away from the Galactic plane in harmonic space as a function of multipole.

• (3)

  We analyse the correlation of dust polarization maps over all four polarized HFI channels from 100 to 353 GHz.

We focus on presenting results of direct relevance to component separation, leaving the astrophysical modelling of the results to follow-up studies. A second paper (Planck Collaboration XII 2020) presents a complementary perspective on dust polarization from an astrophysics perspective, focusing on the statistics of the polarization fraction and angle derived from the 353 GHz Planck maps.

The paper is organized as follows. In Sect. 2, we present the Planck sky maps and their validation. Results from the power-spectrum analysis of the dust polarization maps at 353 GHz are described in Sect. 3. In Sect. 4, the Planck HFI maps in the frequency range 100 to 353 GHz are combined with lower frequency maps from the Planck Low Frequency Instrument (LFI; Planck Collaboration II 2020) and WMAP (Bennett et al. 2013) to characterize polarized foregrounds across microwave frequencies and multipoles, including the correlation between dust and synchrotron polarization. We compare the microwave SEDs of dust polarized emission and total intensity in Sect. 5. We quantify the correlation between Planck-HFI polarized dust maps in Sect. 6. Section 7 summarizes the main results of the paper. The paper also has three appendices. Data simulations used to estimate uncertainties in our data analysis are presented in Appendix A. In Appendix B, we revisit the correlation analysis of the 217- and 353 GHz Planck polarization maps investigated previously in PL, using the PR3 data and E2E simulations. Large tables, are gathered in Appendix C.

2. The Planck PR3 polarization maps

Planck observed the sky in seven frequency bands from 30 to 353 GHz for polarization, and in two additional bands at 545 and 857 GHz for intensity, with an angular resolution from 31′ to 5′ (Planck Collaboration I 2014). The in-flight performance of the two focal-plane instruments, HFI and LFI, are described in Planck HFI Core Team (2011) and Mennella et al. (2011), respectively. For this study, we use the new Planck PR3 maps. The processing of the HFI data is described in Planck Collaboration III (2020) and that of LFI data in Planck Collaboration II (2020).

The 100-, 143-, and 217 GHz HFI maps are made using data from all bolometers, while the 353 GHz maps are constructed using only data from the polarization-sensitive bolometers (PSBs), as recommended in Planck Collaboration III (2020). To characterize the data noise and to compute power spectra at one given frequency that are unbiased by noise, we use maps built from data subsets, specifically the two half-mission and the two odd-even survey maps (Planck Collaboration III 2020)³. In this paper, we focus on results obtained using half-mission maps, but have checked that conclusions would not be changed if we had used odd-even surveys instead. The Planck-HFI data noise and systematics are quantified and discussed in Planck Collaboration III (2020) using the E2E simulations of Planck observations introduced there. The related methodology that we follow to estimate uncertainties from detector noise and residual systematic effects, and to propagate them to the results of our data analysis, is presented in Appendix A.

A posteriori characterization of polarization efficiencies (Planck Collaboration III 2020) suggests small modifications relative to the values used to produce the delivered frequency maps available on the Planck legacy archive. Accordingly, we multiply the PR3 HFI polarization maps at 100, 143, and 217 GHz by 1.005, 0.98, and 1.015, respectively; uncertainties in these factors are of order 0.005. For 353 GHz, no such factor has been determined but we expect it to have the same magnitude as at the other HFI frequencies. Thus, we consider that there is a 1.5% photometric uncertainty on the 353 GHz polarized emission.

In addition, in Sect. 4 we use polarization maps from LFI at 30 GHz, and the K and Ka WMAP channels (Bennett et al. 2013) to separate dust and synchrotron polarized emission and quantify the correlation between the two sources of emission. Because E2E simulations are not available for these data, we compute maps of uncertainties from Gaussian realizations of the data noise. Power spectra of the data noise are derived from the half-difference of half-mission Planck-LFI maps and the difference of year maps for WMAP. We note that it is easy to produce a large number (1000 or more) of data realizations with Gaussian data noise, while only 300 E2E realizations are available for HFI. The data are expressed in thermodynamic (CMB) temperature throughout the paper.

3. Angular power spectra of dust polarization

In this section, we derive angular power spectra of dust polarization from the PR3 maps at 353 GHz. Key improvements in the correction of data systematics allow us to extend earlier work on dust polarization (power spectra and SED) to the lowest multipoles.

3.1. Planck angular power spectra at 353 GHz

The power-spectrum analysis of Planck dust polarization in PXXX was limited to multipoles ℓ >  40, due to residual systematics in the available maps. The improvements made in correcting Planck systematics for the new data release allow us to extend the range of scales over which we can characterize dust polarization.

The EE, BB, TE, TB, and EB power spectra are computed with the XPol code (Tristram et al. 2005). Following the approach in PXXX and PL, to avoid a bias arising from the noise, we compute all of the Planck power spectra using cross-correlations of maps with independent noise, specifically the half-mission maps. To present a characterization of foregrounds that is independent of component-separation methods, we chose not to use the CMB polarization maps described in Planck Collaboration IV (2020). Instead, the CMB contribution is subtracted from the power spectra using the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016). The power spectra shown in the figures and tables below are in terms of , where X ∈ {T, E, B}, Y ∈ {E, B}, and is the XY angular power spectrum. The error bars are derived from the simulations described in Appendix A; they include the cosmic variance of the CMB computed for each sky region, because the CMB is subtracted using the Planck 2015 ΛCDM model.

We examine six nested regions at high Galactic latitude, with an effective sky fraction ranging from 24 to 71%. These regions are defined using the same set of criteria as in PXXX, meant to minimize dust polarization power for a given sky fraction, and with the same apodization (see Fig. 1). The regions differ only in the masking of point sources; we mask a smaller number of sources that are polarized. We keep the same “LRnm” nomenclature, where “nm” is as a percentage. Table C.1 lists other properties of the regions, including the mean specific intensity at 353 GHz, ⟨I₃₅₃⟩ in MJy sr⁻¹, and the mean HI column density, N_H in units of 10²⁰ cm⁻², inferred as in PL from the Planck dust opacity map in Planck Collaboration Int. XLVIII (2016).

Fig. 1. All-sky map showing the sky regions used to measure power spectra, indicated with colours varying from yellow to orange and dark-red. The white region represents the area where the CO line brightness is larger than 0.4 K km s−1, which is excluded from all the sky regions in our analysis. The blue dots represent the areas masked around point sources.

The EE and BB spectra are tabulated in Table C.1 and presented in Fig. 2 for each of our six sky regions. For the lowest multipole bin (ℓ = 2–3), we report a value for only the largest sky region LR71, over which it is best measured. In Fig. 3, we present spectra computed on the northern and southern parts of the LR42, LR52, LR62, and LR71 regions.

Fig. 2. CMB-corrected EE (red diamonds), BB (blue squares), and TE (black circles) power spectra at 353 GHz, for each of the six sky regions that we analyse. The dashed lines represent power-law fits to the data points from ℓ = 40 to 600. The exponents of these fits, αTE, αEE, and αBB, appear on each panel.

3.2. Power-law fits

We performed a χ² fit to the power spectra over the multipole range 40 ≤ ℓ ≤ 600, as in PXXX, using the equation:

(1)

where XY ∈ {EE, BB, TE}. The power-law fits are displayed with dashed-lines in Fig. 2 for the six sky regions and in Fig. 3 for the northern and southern parts of the LR42, LR52, LR62, and LR71 regions. The amplitudes A^EE, expressed in μK², and exponents α_XY are listed in Table 1 for the six sky regions. The exponents are also printed in each panel of Figs. 2 and 3. The error bars on A^EE include a 3% factor from the 1.5% uncertainty on the 353 GHz polarization efficiency.

Fig. 3. Power spectra, as in Fig. 2, but for the northern and southern parts of the LR42, LR52, LR62, and LR71 regions.

The power laws match the fitted data points well, but not perfectly. Indeed, for many regions, including the largest ones with the highest signal-to-noise ratios, the χ² values in Table 1 are larger than the number of degrees of freedom, N_d.o.f. = 24. We note that these χ² values are calculated for exponents fixed at a common value of −2.44. There is evidence for statistically significant variations of the exponents over sky regions. Furthermore, there is a difference between the values for the EE and BB spectra, which for the largest sky region are α_EE = −2.42 ± 0.02 and α_BB = −2.54 ± 0.02, respectively.

Figures 2 and 3 also show the extrapolation of the power laws to low multipoles, which may be compared to the data points at ℓ <  40 not used in the fit. The extrapolation is close to these data points in some cases, but not always.

Dust polarization angular power spectra, like the spectra of synchrotron emission, are related physically to the power spectrum of interstellar magnetic fields. Within the phenomenological models of Ghosh et al. (2017) and Vansyngel et al. (2017), the exponent of the dust power spectrum is found to be close to that of the Gaussian random field used to simulate the turbulent component of the magnetic field. The spectra are expected to flatten towards low multipoles, when the analysis is of an emitting volume sampling physical scales larger than the injection scale of turbulence (Cho & Lazarian 2002). We do not observe such a flattening, but it might well be hidden by systematic variations of the magnetic field orientation over the solar neighbourhood. It will be necessary to extend the work of Vansyngel et al. (2017) to low multipoles in order to assess whether our new results may be accounted for by statistical variance within their model framework.

3.3. Scaling of B-mode power with total intensity of dust emission

In Fig. 4, we plot the amplitude A^BB(ℓ = 80) versus the mean dust total intensity at 353 GHz, ⟨I₃₅₃⟩. The amplitudes are well fit by a power-law of the form ⟨I₃₅₃⟩² (the dashed line in Fig. 4), i.e. with the same exponent as that measured for the amplitude of the total dust intensity angular power spectrum in the far-infrared (Miville-Deschênes et al. 2007) but slightly greater than the value of 1.9 for dust polarization in PXXX. The fit to our results for the six sky regions also matches the measurement reported by Ghosh et al. (2017) for a region of low HI column density in the southern sky with f_sky = 8.5%.

Fig. 4. Scaling of the BB power at ℓ = 80 versus the mean dust total intensity at 353 GHz. The dashed black line is a power-law fit to values for the six sky regions in our analysis (blue dots) with an exponent of two. Also shown are the values for the N–S splits of the regions in Fig. 3 (blue triangles). These results are complemented by the measurement (red diamond) over the southern Galactic cap (fsky = 8.5%) by Ghosh et al. (2017) and that for the BICEP field (black square) after BICEP2 and Keck Array Collaborations (2016).

We also measured the 353 GHz dust B-mode power on the PR3 maps over the BICEP/Keck field using the mask available on the collaboration website. Measurements have been made as well on each of the 300 realizations of the E2E simulations. The dispersion of these measurements provides us with error bars, including both instrumental noise and uncorrected systematics. We find (4.4 ± 3.4) μK² using half-mission maps to compute cross-spectra, and (0.83 ± 3.1) μK² for odd and even surveys. These measurements are consistent with the value derived from the correlation of the Planck 353 GHz PR2 maps with the BICEP/Keck 95- and 150 GHz data in BICEP2 and Keck Array Collaborations (2016); to compare that value, (4.3 ± 1.1) μK² at the reference frequency 353 GHz, with the above results from Planck PSB-only polarization band-integrated data at 353 GHz, we multiply it by the colour correction 1.098² to obtain (5.2 ± 1.3) μK². In the BICEP field ⟨I₃₅₃⟩ = 0.048 MJy sr⁻¹, for which extrapolation of the fit to our measurements gives a signal level of approximately 8 μK². The difference is within the cosmic variance, as estimated by Ghosh et al. (2017) using their statistical model of the dust polarization in the southern Galactic cap.

3.4. Asymmetry between the power in E and B modes

In Table 1, for each of the six regions we list the BB/EE ratio of the amplitudes parameterizing the power-law fits. The weighted mean ratio is BB/EE = 0.524 ± 0.005, a value consistent with that in PXXX. For some regions, but not all of them, we find that the E/B power asymmetry extends to the lowest multipole bins. At low multipoles the measured BB/EE power ratio is in the range of about 0.5–1.

The weighted mean values of the exponents for the EE and BB power spectra are α_EE = −2.39 and α_BB = −2.51, respectively. The weighted dispersions of individual measurements for the six regions are 0.05 and 0.06, respectively. The exponents measured on the northern and southern parts of the LR42, LR52, LR62, and LR71 regions in Fig. 3 fit within this statistical characterization of our results for the full sky regions. The exponents that we find are close to the values reported in PXXX.

However, we find a small difference between the two exponents, which suggests that the asymmetry changes slightly as a function of multipole. Such a difference is not unexpected. The filamentary structures in the cold neutral interstellar medium have mainly E-mode polarization, due to the statistical alignment of the magnetic field orientation with matter (Clark et al. 2014; Planck Collaboration Int. XXXVIII 2016; Ghosh et al. 2017).

3.5. The TE correlation

Planck Collaboration Int. XXXVIII (2016) related the TE correlation to the observed alignment between filamentary structures and the magnetic field in the diffuse ISM, while Caldwell et al. (2017) discussed it theoretically in the context of MHD turbulence. However, the new data shown here in Figs. 2 and 3 show that the TE correlation extends down to the lowest multipoles, which characterize dust polarization on angular scales larger than those of interstellar filaments. To examine this further, we performed χ² fits of a power law to the TE spectra, as for EE and BB, over the multipole range ℓ = 40 − 600. The parameters of the fits are listed in Table 1 and displayed in Figs. 2 and 3. The data points at ℓ <  40, not included in the fit, are close to the extrapolation of the power laws to low multipoles.

These new results show that the filamentary structure of the magnetized interstellar medium alone cannot account for the observed TE correlation. At least for the lowest multipoles, the correlation must have another origin that will need to be explored in future studies. One possibility is that the low-ℓ TE correlation arises from the correlation between the local structure of the GMF with the geometry of the Local Bubble cavity (Alves et al. 2018).

The weighted mean value of the exponent is α_TE = −2.49, slightly different than α_EE = −2.39. The TE spectrum is shallower (i.e. the absolute value of α_TE is smaller) than that measured on average for HI column density maps (Miville-Deschênes et al. 2003; Martin et al. 2015; Blagrave et al. 2017). However, using line profile decomposition to isolate gas with the lowest velocity dispersion (the cold neutral medium or CNM), Martin et al. (2015) and Ghosh et al. (2017) provide evidence that the angular power spectrum of the column density of the CNM gas is shallower, in particular with exponent about −2.4 in the extended SGC34 region defined by the latter (a 3500 deg² region comprising 34% of the southern Galactic cap with ). As quantified by the modeling in Ghosh et al. (2017), this is in agreement with the idea that the TE correlation, and the E/B asymmetry, at ℓ >  40 are related to the statistical alignment of the magnetic field with filamentary structure in the cold medium (Clark et al. 2015; Planck Collaboration Int. XXXVIII 2016; Kalberla et al. 2016).

Table 1 gives values of the ratio of the amplitudes of the TE and EE power spectra. The weighted mean value of the TE/EE ratios is 2.76 ± 0.05. We also combine the dust TE, EE, and TT spectra at 353 GHz to compute the dimensionless correlation ratio discussed by Caldwell et al. (2017) and introduced in the context of the CMB in Appendix E.3 of Planck Collaboration XI (2016). The ratio is plotted versus multipole in Fig. 5 for the six regions. The weighted mean of all measurements for all sky regions and multipole bins is . The data show significant scatter, but no systematic dependence on multipole down to the lowest ℓ bins or on the sky region.

Fig. 5. TE correlation ratio versus multipole. The data points are plotted using distinct symbols and colours (see legend at the top) for each of the six sky regions. The error bars are derived from the E2E simulations.

3.6. TB and EB power spectra

The TB and EB angular power spectra are presented in Fig. 6. We find a positive TB signal. A similar result was reported using earlier Planck data in PXXX. On the largest sky regions providing the best signal-to-noise ratio, the power ratio TB/TE is about 0.1 from a power-law fit (exponent fixed at −2.44) over the ℓ = 40–600 multipole range. The correlation ratio , about 0.05, is also much lower than . The EB signal is consistent with zero. The EB/EE power ratio is smaller than about 0.03.

Fig. 6. Power spectra of TB (red diamonds) and EB (blue squares) at 353 GHz for the six sky regions. The error bars are derived from the E2E simulations. A power-law fit to the TB data (solid red line) reveals an overall positive TB signal, not seen in the E2E simulations. The EB power (solid blue line fit) is consistent with zero.

The E2E simulations in this paper allow us to check that the TB power does not arise from a known systematic error. For example, a systematic error in the orientation of the Planck bolometers at 353 GHz would induce leakage of the TE power to TB and from the EE and BB power to EB (Abitbol et al. 2016). To account for a ratio TB/TE = 0.1, the error would need to be 3°, a value that is one order of magnitude larger than the uncertainties on the orientation of the HFI PSBs determined from CMB data analysis for the 100, 143, and 217 GHz channels (see Appendix A.6 in Planck Collaboration Int. XLVI 2016).

We do not see any systematic effects that could produce the TB signal. If it is indeed real, this indicates that the dust polarization maps do not satisfy parity invariance. Although there is no reason for Galactic emission to preserve mirror symmetry, to our knowledge there is no straightforward interpretation of this observed asymmetry. The TB signal, at low multipoles, might arise from the structure of the mean magnetic field in the solar neighborhood. It might also be related to reference quantities of magnetohydrodynamic turbulence that are not parity invariant, such as the magnetic helicity (the volume integral of the scalar product between the vector potential and the magnetic field; see e.g. Blackman 2015) and/or the cross-helicity (the integral of the scalar product between the gas velocity and the magnetic field; see e.g. Yokoi 2013). These possible links will need to be explored in further studies.

Within the context of CMB experiments, as discussed in Abitbol et al. (2016) a non-zero dust TB signal can limit the accuracy to which TB and EB spectra at microwave frequencies may be used to check the orientation of the polarimeter.

4. Dust and synchrotron polarized emission at microwave frequencies

We now calculate cross-power spectra, build models for them, and compare the foreground signals to the CMB. Specifically, in Sect. 4.1 using cross-spectra we characterize Galactic polarized emission, including the correlation between dust and synchrotron polarization, as a function of frequency and multipole. In Sect. 4.2, we fit these data with a spectral model and present the parameters determined. Galactic polarized foregrounds as quantified here are compared to the CMB primordial E- and B-mode signals as a function of frequency and multipole in Sect. 4.3.

4.1. Cross-power spectra

For this study, we consider single and inter-frequency cross-spectra among the four polarized channels of Planck-HFI, at 100, 143, 217, and 353 GHz, as well as the lowest frequency channel of Planck-LFI at 30 GHz, and the two lowest frequencies of WMAP at 23 and 33 GHz. The three channels of LFI and WMAP provide the highest signal-to-noise ratio on synchrotron polarization; we use them to estimate the synchrotron contribution to the lowest HFI frequencies and characterize the spatial correlation between polarized dust and synchrotron sources of emission.

Single frequency cross-spectra are computed using maps with independent statistical noise made with data subsets, to avoid noise bias. For Planck-HFI, we use the half-mission maps. For Planck-LFI, we separate data from even and odd years. For WMAP, we combine the first four years on the one hand and the subsequent five years on the other hand. For inter-frequency cross-spectra, we consider all the possible combinations among the frequency channels being used. In total, we obtain 21 cross-spectra that combine observations at two distinct frequencies and 7 cross-spectra at a single frequency. The uncertainties on power spectra are again computed from E2E simulations, as described in Appendix A.

All 28 spectra are computed for each of the six sky regions described in Sect. 3.1, within nine multipole bins in the range 4 ≤ ℓ ≤ 159. The specific multipole bins are top-hat (flat) in the following ranges: 4–11; 12–19; 20–39; 40–59; 60–79; 80–99; 100–119; 120–139; and 140–159. Low signal-to-noise ratios prevent us from deriving meaningful SED parameters at higher multipoles. Figure 7 presents an example for B modes in the LR62 region for two multipole bins, ℓ = 4–11 and 40–59.

Fig. 7. BB cross-spectra 𝒟ℓ(ν1 × ν2) versus the effective frequency νeff = (ν1 × ν2)0.5, for the LR62 sky region and two multipole bins: ℓ = 4-11 (top plot) and 40–59 (bottom). Yellow and blue colours represent data values from single and inter-frequency cross spectra, respectively. Bottom panel: within each plot shows the residuals from the fits normalized to the 1 σ uncertainty of each data point. Lower frequency data (left) points are dominated by the SED of synchrotron polarized emission, while higher frequency (right) data characterize dust polarized emission, and those at the centre characterize the correlation between the two sources of emission. Differences between the two plots illustrate that both the ratio between synchrotron and dust power and the correlation between these two sources of polarized emission decrease for increasing multipoles.

4.2. Spectral model

Our SED analysis includes polarized synchrotron emission spatially correlated with polarized thermal dust emission (Kogut et al. 2007; Page et al. 2007; Planck Collaboration Int. XXII 2015; Planck Collaboration X 2016). We use the following spectral model, introduced by Choi & Page (2015):

(2)

where X ∈ {E, B} and is the amplitude of the XX cross-spectrum between frequencies ν₁ and ν₂ (expressed in GHz) within a given multipole bin ℓ, expressed in terms of brightness temperature squared. The Planck function B_ν(T_d) is computed for a fixed dust temperature T_d = 19.6 K, derived from the fit of the SED of dust total intensity at high Galactic latitude in Planck Collaboration Int. XXII (2015). We use a fixed temperature because, over microwave frequencies, the dust SED depends mainly on the dust spectral index of the modified black-body (or MBB) emission law and the temperature cannot be determined independently of the spectral index. As discussed in Planck Collaboration Int. XXII (2015) and Choi & Page (2015), the cross-correlation between dust and synchrotron polarization might arise from the magnetic field structure but might also include a contribution from variations of the synchrotron spectral index and anomalous microwave emission (AME) if it is polarized (Hoang & Lazarian 2016a; Draine & Hensley 2016; Génova-Santos et al. 2017).

The spectral model has five parameters: the two amplitudes A_s and A_d and the two spectral indices β_s and β_d, characterizing the synchrotron and dust SEDs, respectively; and the correlation factor ρ quantifying the spatial correlation between synchrotron and dust polarized emission. In Eq. (2), the synchrotron and MBB emission are expressed in terms of brightness temperature, whereas the data are in thermodynamic units. The conversion between the two is accomplished by two factors. The first, U, is a unit conversion from the thermodynamic units to brightness temperature units for some adopted reference spectral dependence, performing the appropriate integrations over the bandpass. The second, C, is a colour correction from the actual spectrum of the model to the adopted reference spectral dependence, again with bandpass integrations. Accordingly, the spectrum is converted into units of the data by multiplication by C/U, and in the application to the fit of the spectral model in Eq. (2) by multiplication by (C/U)₁(C/U)₂. These factors were computed as in Planck Collaboration Int. XXII (2015), for Planck using the procedures hfi_unit_conversion and hfi_colour_correction (for both HFI and LFI) and the instrument data files described in the Planck Explanatory Supplement⁴, and for WMAP the formulae and tabulations in Jarosik et al. (2003). Here, for both HFI and LFI the adopted reference spectral dependence is I_ν ∝ ν⁻¹ (see discussion in Planck Collaboration IX 2014 and the Planck Explanatory Supplement⁵), whereas for WMAP it is constant Rayleigh-Jeans temperature. By construction, the ratio C/U does not depend on the adopted choice. The conversion factors used are listed in Table 2. These are very close to the factors in Table 3 of Planck Collaboration Int. XXII (2015), though here at 353 GHz the evaluation is for the PSBs only. The values of C are evaluated for the following SED. For the LFI and WMAP channels used, the synchrotron component dominates, for which we assume β_s = −3, while for the Planck HFI channels the polarized dust MBB spectrum dominates, for which we assume β_d = 1.5 and T_d = 19.6 K.

We fit our spectral model to the EE and BB spectra separately, for each sky region and for each multipole bin independently. Before fitting, we subtract the amplitude of the CMB power spectrum, estimated from the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016), from each data point. The fit is carried out in two steps. First, we fit the model of Eq. (2) using the MPFIT code, which uses the Levenberg-Marquardt algorithm to perform a least-squares fit. We then compute the weighted mean and standard deviation of β_s over the MPFIT results for all sky regions and multipole bins, finding β_s = −3.13 ± 0.13. This value of β_s is consistent with those obtained by Fuskeland et al. (2014) and Choi & Page (2015) using all frequency channels of WMAP, as well as that, −3.22 ± 0.08, reported by Krachmalnicoff et al. (2018) for the frequency range 2.3 − 33 GHz, combining data from the S-band Polarization All-Sky Survey (S-PASS) at 2.3 GHz, WMAP, and Planck. We use it as a Gaussian prior for a second fit of the same data with the same model. This second fit is performed with a Monte Carlo Markov chain (MCMC) algorithm. In both fits we assume that the data points are independent. We checked on the E2E realizations that this is true for the B-mode data. For E-mode, the CMB variance introduces a slight correlation that we neglect. We adopt this two-step procedure because when attempting to fit β_d without a prior on β_s we found spurious results for a few combinations of ℓ_bin and sky regions, when the signal-to-noise ratio in the low-frequency channels is too low to constrain the synchrotron SED adequately.

An example is given in Fig. 7, also showing the residuals from the fit. The χ² values for all fits are listed in Tables C.2 and C.3 for the EE and BB spectra, respectively. The results obtained on the simulated maps (Fig. A.4) show that the fit parameters match the input values without any bias.

Continuing the example, Fig. 8 shows the posterior distribution of the model parameters obtained through the MCMC algorithm, for BB data, the LR62 region, and the ℓ = 40–59 bin. Best-fit parameters are computed as the median value of the posterior distributions, while errors are obtained from the 16th and 84th percentiles (68% confidence interval). For all regions and multipole ranges, values for A_d, A_s, β_d, β_s, and ρ are listed in Tables C.2 (EE) and C.3 (BB).

Fig. 8. Posterior distribution for each of the parameters of the spectral model in Eq. (2), as obtained through the MCMC fitting algorithm for BB data points. The MCMC results illustrated here are for the LR62 region and the multipole bin ℓ = 40–59, one of the two cases shown in Fig. 7. The diagonal shows the probability distribution of each parameter. Median values are As = 0.6 ± 0.1, Ad = 137 ± 2, βs = −3.15 ± 0.17, βd = 1.50 ± 0.02, and ρ = 0.17 ± 0.04.

We do not list the amplitudes A_d and A_s of the dust and synchrotron emission but we note that as expected values of A_d are close to the values of the amplitudes in Table C.1. In Fig. 9, A_d and A_s for EE and BB are plotted versus multipole for the six sky regions. As in the spectra for each region in Fig. 2, A_d has a power-law dependence on ℓ and a systematic increase with (see e.g. Fig 4) that applies down to lower multipoles beyond ℓ = 40. On the other hand, for the multipole bin 4–11 the B-mode synchrotron amplitude is roughly constant over the six sky regions. As a corollary, for this multipole bin the ratio between dust and synchrotron B-mode polarization increases by about one order of magnitude from the smallest sky region, LR24, to the largest one, LR71. We point out that this result is specific to our set of sky regions, which are defined using the dust total intensity map to minimize dust power for a given sky fraction.

Fig. 9. Amplitudes of EE and BB power spectra for dust and synchrotron emission at 353 and 30 GHz, respectively, shown for each sky region and each multipole bin. The As and Ad parameters of our spectral model from Eq. (2) are thermodynamic (CMB) temperature in μK2. Where the synchrotron amplitude is compatible with zero at the 1 σ level, we report an upper limit on As (68% confidence limit) with triangles pointing down.

Krachmalnicoff et al. (2018) have characterized the synchrotron polarized foreground emission analysing maps of the southern sky from S-PASS at 2.3 GHz. Comparison with our synchrotron results in Fig. 9 is not immediate because power spectra are not measured over the same sky regions. Further, the signal to noise ratio of the S-PASS data for synchrotron emission is larger than that of WMAP and Planck, which is a critical advantage in characterizing the faint polarization signal at high Galactic latitude. However, contamination by Faraday rotation is likely to be significant for their largest sky regions extending down to Galactic latitude |b| = 20°.

Figure 10 plots the two parameters ρ and β_d (not β_s because of the prior applied) for EE and BB. The top panels show that ρ, which quantifies the correlation between dust and synchrotron polarization, decreases with increasing multipole and is detected with high confidence only for ℓ ≲ 40. The correlation might extend to higher multipoles, but the decreasing signal-to-noise ratio of the synchrotron polarized emission precludes detecting it. These results are consistent with the analysis done by Choi & Page (2015) using all frequency channels of WMAP. The bottom panels show that the spectral index β_d has no systematic dependence on multipole or sky region, except for the lowest multipole bin. The dust spectral indices are further discussed in Sect. 5.

Fig. 10. Fit parameters ρ and βd for E- and B-mode polarization versus multipole. Open symbols for ρ represent the cases where the synchrotron amplitude is compatible with zero, making it difficult to measure the correlation.

4.3. Foregrounds versus CMB polarization

Next, Galactic foregrounds are compared to CMB E- and B-mode polarization to quantify the challenge of component separation for measuring the low-multipole E-mode CMB signal from reionization (Fig. 11), and also for detecting primordial B modes (Figs. 12 and 13). The results of our spectral analysis allow us to update earlier studies (see e.g. Dunkley et al. 2009; Krachmalnicoff et al. 2016; Planck Collaboration X 2016).

Fig. 11. Dust and synchrotron E-mode power versus multipole. The dust power at 95 and 150 GHz and that of synchrotron at 95 GHz are compared with the CMB E-mode signal (red-line) computed for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016) and a Thompson scattering optical depth τ = 0.055 from Planck Collaboration Int. XLVI (2016). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018).

Fig. 12. Dust and synchrotron B-mode power versus multipole. The dust power at 95 and 150 GHz, and that of synchrotron at 95 GHz are compared with CMB B modes from primordial gravitational waves (grey lines) for three values of the tensor-to-scalar ratio, r = 0.1, 0.01, and 0.001, and from lensing (blue line) for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018).

Fig. 13. Dust and synchrotron B-mode power versus frequency for two multipole bins: ℓ = 4–11 (top) and 60–79 (bottom). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018). The primordial CMB B-mode signal, averaged within the appropriate ℓ bin, is plotted with dashed lines for three values of the tensor-to-scalar ratio: r = 0.1; 10−2; and 10−3. The solid line represents the lensing B-mode signal for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016).

To prepare Figs. 11 and 12, we use the results of our spectral fitting to compute the dust and synchrotron E- and B-mode power at frequencies 95 and 150 GHz, which correspond to the two microwave atmospheric windows providing the best signal-to-noise on the CMB for ground-based observations. In both figures, the dust power is represented by a coloured band that spans the signal range from the smallest (LR24) to the largest (LR71) sky regions in our analysis; the lower and upper edges of the band represent power-law fits of the values of A_d listed in Tables C.2 and C.3. For synchrotron and LR71, we apply the same procedure fitting A_s values. For LR24, the signal to noise ratio of our synchrotron results is too low to compute a reliable fit. We choose instead to plot the results from the S-PASS data in Krachmalnicoff et al. (2018) for their smallest (|b| > 50°) sky region. The scaling of the power spectrum amplitude from 2.3 to 95 GHz is done using their determination of the spectral index (β_s = −3.22).

The dust E-mode power at 95 and 150 GHz and that of synchrotron at 95 GHz are compared with the CMB, as a function of multipole, in Fig. 11. Similarly, in Fig. 12, the B-mode foregrounds at the same two frequencies are compared with the CMB primordial and lensing signals. The primordial B-mode signal has two broad peaks in two multipole ranges, ℓ = 2–8 and 30–200, associated with reionization and recombination, respectively, the amplitude of which scales linearly with the tensor-to-scalar ratio r. The E- and B-mode reionization bumps at low multipoles are computed here for a Thompson scattering optical depth τ = 0.055 from Planck Collaboration Int. XLVI (2016).

Figure 12 shows that the synchrotron power decreases more steeply than the dust power with increasing ℓ. Consequently, polarized synchrotron is a more significant foreground for the reionization peak than for the recombination peak.

In Fig. 13, the dust and synchrotron BB power is plotted versus frequency for two multipole bins ℓ = 4–11 (top plot) and 60–79 (bottom plot), which roughly correspond to the reionization and recombination peaks of the primordial B-mode CMB signal, respectively. The lower and upper edges of the dust band are drawn combining A_d values with spectral indices β_d, both listed in Table C.3, for the LR24 and LR71 sky regions. For synchrotron, as in Figs. 11 and 12, we use the results from the S-PASS data in Krachmalnicoff et al. (2018) for their smallest (|b| > 50°) sky region to draw the lower edge of the coloured band. The two polarized foregrounds have comparable amplitudes at a frequency that depends on the multipole and the sky region. For average ℓ_bin = 7.5 (top plot) the amplitudes are equal at ∼75 GHz for both the lower and upper edges of the bands, whereas for ℓ_bin = 69.5 (bottom) this equality occurs at a lower frequency ∼60 GHz. For higher frequencies, dust quickly dominates synchrotron. For example, for ℓ_bin = 69.5, the BB dust and synchrotron signals are equal at 60 GHz, while at 100 GHz the dust and synchrotron powers differ by two orders of magnitude.

Our analysis stresses the accuracy with which dust and CMB B modes must be separated to search confidently for primordial B modes down to r = 0.01. At this sensitivity level for sub-orbital experiments targeting the recombination peak at 95 and 150 GHz, e.g. the BICEP/Keck Array ground-based experiment (BICEP2 and Keck Array Collaborations 2016) and the Spider balloon-borne experiment (Fraisse et al. 2013), synchrotron polarization appears not to be a significant foreground over the relevant high latitude southern sky areas at |b| > 50° used to draw the lower edge of the band. However, the exact level of contamination will depend in detail on the sky region observed and how the synchrotron power extrapolates from 2.3 GHz there.

5. Microwave SED of polarized dust emission

This section focusses on the microwave SED of dust emission that is of interest for component separation and as a constraint on dust emission models.

5.1. Spectral index of dust polarized emission

Within the approximation of an MBB emission law and given a dust temperature, the microwave SED of dust emission is determined by the value of the dust spectral index, β_d. This index parameterizes the separation of the dust and CMB components and the Planck data constrain it better than ground-based data thanks to Planck’s 353-GHz channel.

We compute the mean values and for E- and B-mode polarization from the results of the spectral fitting from Sect. 4 in Tables C.2 and C.3. The uncertainty-weighted average of the differences between and , computed over all multipole bins and sky regions, is . We consider the significance of this difference to be marginal because the statistical error-bar assumes that the measurements for the different sky regions are independent. Averaging differences for the LR71 region alone, we find .

The difference between and is small and so we averaged them. Specifically, the uncertainty-weighted average of the fit results for all multipole bins and sky regions is , where the error bar includes the uncertainty from the polarization efficiencies of HFI (Sect. 2) and the uncertainty from the CMB subtraction, which affects the determination of . This is the uncertainty of the mean; the weighted dispersions of individual measurements are 0.046 and 0.034 for E and B modes, respectively. This value of is lower than the mean polarization index 1.59 ± 0.02 derived from the analysis of earlier (PR2) Planck data (Planck Collaboration Int. XXII 2015). This difference reflects correction of data systematics between the PR2 and PR3 polarization maps. We checked that it does not come from the data analysis by analysing the PR2 data in the same way as the PR3 data in this paper.

5.2. Dust polarization SED from blind component separation

The dust SEDs for E- and B-mode polarization were determined jointly with the corresponding synchrotron SEDs using the SMICA (Spectral Matching Independent Component Analysis) method of blind component separation described in Cardoso et al. (2008), Planck Collaboration IX (2016) and Planck Collaboration IV (2020). In brief, the method consists of fitting all of the auto- and cross-spectra from 30 GHz to 353 GHz to a model consisting of a superposition of the CMB, two foreground emission components, and noise. The fit is performed under very mild constraints, the free parameters being the angular spectrum of the CMB, the SED of each foreground emission component (assumed independent of angular scale), the angular spectra of each foreground emission component and their cross-spectrum, and the noise spectrum at each frequency. No prior spectral models of the SEDs are assumed; we do not assume that the dust SED is an MBB or that the synchrotron SED is a power law.

Fitting such a model determines, at the spectral level, a unique global foreground contribution that corresponds to two underlying templates. However, because the model allows for an arbitrary angular correlation between those two templates, as well as an arbitrary SED for each of them, the templates are linearly degenerate, meaning that each can be an arbitrary linear combination of synchrotron and dust emission. We choose to resolve this degeneracy by selecting the (essentially unique) linear combinations, such that one template has no contribution at 353 GHz while the other has no contribution at 30 GHz. The latter corresponds to the dust foreground.

The SMICA component separation was performed over the LR71 sky region for comparison with our data analysis. The resulting dust SEDs for E- and B-mode polarization are presented in Fig. 14. These SEDs, coming from blind component separation, are remarkably close to a single-temperature MBB over the full range of Planck polarization observations, despite the fact that an MBB spectral shape was not a prior assumption. Performing MBB fits after the fact to the SMICA dust spectral data in Fig. 14 (again with T_d = 19.6 K and using colour corrections as described in Sect. 4.2), we find a mean spectral index of  = 1.53 ± 0.02, taking into account the 1.5% uncertainty on the polarization efficiency at 353 GHz. The E- and B-mode data intensities, each normalized to 1 at 353 GHz, and uncertainties are listed in Table 3. For comparison, we also list the corresponding values for a MBB SED with  = 1.53 ± 0.02.

Fig. 14. Dust SEDs for E- and B-mode polarization derived from the SMICA component-separation procedure (Planck Collaboration IV 2020). The two grey lines represent MBB fits to the E- (red diamonds) and B-mode (blue squares) data points with a temperature of 19.6 K. The polarization spectral index derived from the fits is  = 1.53 ± 0.02. The residuals to each fit, normalized to the 1 σ data uncertainty, are plotted in the lower panel.

The fit is in excellent agreement with our determination in Sect. 5.1. This agreement is perhaps not that surprising because our approach to the data analysis is in some aspects quite similar to that used by SMICA. In both cases, the foreground SEDs are determined by fitting cross-spectra. Both methods allow for correlation between the two foreground components. However, the two methods differ in their simplifying assumptions. We constrain the dust and synchrotron SEDs to be the MBB and power-law parametric models, while SMICA assumes that the SEDs are scale invariant. The agreement of the SEDs is reassuring and a cross-validation of the assumptions, as well as of the technical implementation.

The BB/EE power ratio from SMICA is 0.60, whereas we find BB/EE = 0.53 ± 0.01 (Table 1). The slightly higher BB/EE power ratio could result from the fact that the BB/EE power ratios in our analysis are determined at ℓ ≥ 40, while SMICA includes lower multipoles. When further constrained to a multipole range approximating ours, the ratio is 0.57.

5.3. Difference between spectral indices for polarization and total intensity

The spectral model in Eq. (2) cannot be applied to the TT spectra because in addition to synchrotron and dust thermal emission there are two other Galactic components, namely AME and free-free emission, that contribute to the total intensity of the Galactic signal (Gold et al. 2011; Planck Collaboration X 2016; Planck Collaboration XXV 2016). To compare the SEDs of dust polarization and total intensity, we follow a method similar to that used in Planck Collaboration Int. XXII (2015) correlating emission in the 217- and 353 GHz HFI channels. We work in harmonic space to assess any SED dependence on multipole and to be able to compare these results to those from the SED fitting. In doing this, we implicitly assume that AME and free-free may be neglected at these two frequencies.

We compute the colour ratio,

(3)

for the TT, EE, and BB spectra. The ratios are colour corrected, as described in Sect. 4.2. We derive the corresponding spectral indices for a dust temperature of 19.6 K. To compute , we subtract CMB anisotropies using the map produced with the SMICA component-separation method. The 353 GHz power spectra are computed using half-mission data subsets.

The spectral indices are listed for each sky region and multipole bin in Table C.4 for the Planck PR3 data. The results are also presented in Fig. 15. The sky emission model that we use for simulating the total intensity maps includes anisotropies of the cosmic infrared background (Planck Collaboration XXX 2014). For the simulations, we retrieve the dust spectral indices adopted as input (1.50 for the total intensity and 1.59 for polarized intensity) with no bias.

Fig. 15. Comparison of spectral indices of dust polarized emission and total intensity. The spectral indices are derived from the 353-to-217 GHz colour ratio. Plots to the left show the results obtained from our simulated maps, and the ones to the right are from the Planck data. Distinct symbols are used to represent each of the six sky regions, as in Fig. 10. For the simulations, the dashed lines represent the input dust spectral indices (, ). For the data, the dashed lines represent the mean measured dust spectral indices (, ).

For the Planck maps, the dust spectral index for polarized intensity averaged over all regions and all ℓ bins is , taking into account the 1.5% uncertainty on the polarization efficiency at 353 GHz. This value agrees well with that inferred from the multi-frequency spectral analysis in Sects. 5.1 and 5.2 above. The corresponding value for total intensity is , with much smaller uncertainty⁶. The spectral indices for polarization and total intensity differ by 0.05 ± 0.03. This difference is smaller than that reported in Planck Collaboration Int. XXII (2015) analysing earlier Planck data.

We checked the consistency of our derivation of the dust spectral index for polarization with the component separation methods in Planck Collaboration IV (2020), by computing maps of from Planck 217 and 353 GHz CMB-subtracted maps smoothed to a 3° beam. This is illustrated in Fig. 16, where the probability distribution of the 217–353 GHz colour ratio for dust polarized intensity, computed over the LR71 sky region, is shown for each of the component separation methods in Planck Collaboration IV (2020). For all methods, the median value of , inferred from the colour-ratio for a dust temperature of 19.6 K, is consistent with our estimate 1.53 ± 0.02 in Sects. 5.1 and 5.2. We point out that the scatter in measured colour-ratios is dominated by data noise.

Fig. 16. Illustration of the consistency between our analysis and component separation methods. The probability distribution of the 217–353 GHz colour ratio for dust polarized intensity, computed over the LR71 sky region from Planck CMB-subtracted maps smoothed to a 3° beam, is plotted for each of the component separation methods in Planck Collaboration IV (2020). The vertical line is the value derived from our analysis. For the unit conversion factors and color corrections, and our modification of the 217 GHz polarization efficiency, it corresponds to the spectral index from Sects. 5.1 and 5.2. The width of the line represents the error bar.

5.4. Impact on dust modelling

These results from the spectral fitting of the polarized dust SED provide an additional constraint for dust modelling. Reviewing the spectral fit in Sect. 4, for ℓ ≤ 100 all of the χ² values of the spectral fit (listed in Tables C.2 and C.3) are lower than the number of degrees of freedom. Therefore, to the sensitivity of the Planck data, a single temperature MBB emission law is a satisfactory model of the polarized dust emission. This same conclusion is supported by the further analyses in the subsections above. There is no evidence for a flattening or steepening of the dust SED, which could in principle result from a variation of spectral index with frequency as reported from laboratory studies of silicate grains (Demyk et al. 2017), or from a significant contribution from magnetic dipolar emission from magnetic nano-particles (Draine & Hensley 2013).

Interstellar dust is often modeled as a mixture of silicates and carbon grains (e.g. Li & Draine 2001; Draine & Fraisse 2009; Compiègne et al. 2011; Jones et al. 2013; Siebenmorgen et al. 2014; Guillet et al. 2018). A difference between and might be evidence that these two dust components have distinct spectral indices and polarization properties. However, the difference that we have found is small and not of high statistical significance. This result suggests that the emission from a single grain type dominates the long-wavelength emission in both polarization and total intensity. If the emission from silicate grains dominates that of carbon grains in polarization – as it is often assumed (Andersson et al. 2015) – this should also hold for the total dust intensity at long-wavelengths.

The alignment of interstellar silicates may be effective irrespectively of whether the grains contain magnetic inclusions (Lazarian & Hoang 2007; Andersson et al. 2015; Hoang & Lazarian 2016a). If silicates do have magnetic inclusions, or if interstellar dust comprises free-flying magnetic grains, the microwave dust emission may include a significant contribution from magnetic dipole emission (Draine & Lazarian 1999; Draine & Hensley 2013; Hoang & Lazarian 2016b; Hensley & Bull 2018). The close match between and constrains this contribution. More generally, the dust polarization SEDs in Table 3 may be used in combination with the dust total intensity SED in Planck Collaboration Int. XXII (2015) (corrected for the 1.5% upward photometric calibration change from the PR2 to PR3 data at 353 GHz) to test dust emission models. This detailed comparison between data and models is beyond the scope of this paper.

6. Correlation of dust polarized emission across microwave frequencies

Interstellar processes couple the emission properties of dust and grain alignment with the density structure of matter and that of magnetic fields (Hoang & Lazarian 2016a; Fanciullo et al. 2017). Likewise, the cosmic-ray energy spectrum, and thereby the synchrotron emission spectrum, depend on the magnetic field structure (Strong et al. 2011). These physical couplings break the simplest assumption for component separation, by which the spectral frequency dependence of the Galactic polarization and its angular structure on the sky are separable (Tassis & Pavlidou 2015; Poh & Dodelson 2017). The couplings make polarized foregrounds intrinsically complex, in ways that have yet to be characterized statistically for optimizing the component separation and taking into account Galactic residuals in the CMB likelihood function. This is a critical issue for the analysis of CMB polarization because spatial variations of the spectral behaviour of polarized dust emission can mistakenly be interpreted as a (false) detection of primordial CMB B modes.

PL analysed the correlation between the HFI dust polarization maps at 217 and 353 GHz. In Appendix B, using the new Planck maps, we update and extend the PL analysis (Appendix B.1). Uncorrected systematics and correlated noise in the data limit how tightly the decorrelation can be constrained. However, these effects change with frequency and so can potentially be mitigated by analysis across many frequencies. In Sect. 6.1, we present such a multi-frequency correlation analysis, making use of the four polarized HFI channels from 100 to 353 GHz. The implications of this new analysis of the Planck data for on-going and future CMB B-mode experiments are discussed in Sect. 6.2.

6.1. Multi-frequency correlation analysis of dust polarization

The spectral model introduced in Sect. 4.2 assumes that the dust and synchrotron polarized emission signals are each perfectly correlated across microwave frequencies. To test this hypothesis, we repeat the spectral fitting with a model modified to allow for a loss of correlation for dust polarization. The dust contribution to the amplitude of BB cross-spectra between frequencies ν₁ and ν₂ is

(4)

where the frequencies ν₁ and ν₂ are expressed in GHz and the adopted function f_d from Appendix B of Vansyngel et al. (2017) is

(5)

The loss of correlation introduced by the parameter δ_d increases with the frequency ratio ν₁/ν₂. From δ_d we also re-express the decorrelation in terms of the spectral correlation ratio (see Eq. (B.1)) for comparison with the two-frequency results presented in PL and in Sheehy & Slosar (2018), and for the PR3 data in Appendix B.

We fit this model over the four HFI polarized Planck frequencies 100, 143, 217, and 353 GHz, for the six sky regions LR24 to LR71. Synchrotron polarization is ignored because it is negligible in this frequency range (Sect. 4.3). We carry out this analysis for the BB cross-spectra computed from the Planck data and the E2E simulations, for the multipole range ℓ = 50–160, relevant to the search for primordial B modes at the recombination peak. To allow readers to fit an alternative spectral model, we list the data values and uncertainties for the LR71 sky region in Table 4. We also provide the corresponding values for the 300 E2E simulations as a FITS Table, which may be used to assess the significance of such an alternative analysis.

We perform an MCMC fit to the Planck data and to the mean of the E2E simulations computed over the 300 E2E realizations. The uncertainties are in both cases inferred from the dispersion of spectra computed with the E2E simulations. In Fig. 17, we show for the LR62 region the posterior probability distribution of the model parameters A_d, β_d, and the correlation ratio inferred from δ_d. The values of the model parameters are listed in Table 5 for the data and the mean of the simulations for all six regions. The dust sky model used in the simulations has a perfect correlation across frequencies (Appendix A.2), that is for this dust model, δ_d = 0 and . The values of in Table 5, inferred from the best-fit value of δ_d for the mean of the 300 E2E realizations, are consistent with 1 within a fraction of the 1 σ error bars, for all sky regions. This result shows that there is no bias introduced by neglecting the synchrotron contribution at 100 GHz, even though it is present in the FFP10 sky model (Appendix A.1). In this model, the contribution of synchrotron to the BB power at 100 GHz, in the multipole bin ℓ = 50–160, rises from 4 to 19% for decreasing f_sky from LR71 to LR24.

Fig. 17. Posterior distribution for each of the parameters of the spectral model with decorrelation given in Eq. (5), as obtained through the MCMC fitting algorithm for BB data points. The MCMC results illustrated here are for the LR62 region and the multipole range 50–160. Median values are Ad = 97.1 ± 1.2, βd = 1.54 ± 0.02, and .

We obtain histograms of parameter values, fitting the spectral model in Eq. (5) to each of the 300 E2E realizations. To do this, we use the least-squares MPFIT algorithm because the MCMC fit is too computationally-intensive to be run 300 times. We checked that the two methods provide consistent parameter values for the Planck data and for the mean of the E2E simulations. The probability distributions of inferred from δ_d values measured on the E2E realizations for each sky region are presented in Fig. 18. Lower limits on from the E2E simulations are listed in Table 5. These are based on the 95% confidence interval, thus on the 2.5th percentile of the histograms.

Fig. 18. Distribution of the correlation ratios inferred from δd on the six sky regions for the ℓ range 50–160. The histograms are computed from the 300 E2E simulations using half-mission data splits. The dashed lines represent the median values on each sky region. This median value, μ, and the standard deviation, σ, are printed in the upper right of each panel. The lower limits on in Table 5 are derived from the 2.5th percentile of the distribution for each sky region.

The limits from the multi-frequency analysis are tighter than the corresponding ones in Table B.1, derived from the 217- and 353 GHz correlation alone (see Appendix B and for convenience reproduced in Table 5). However, it is important to keep in mind that the limits derived from our multi-frequency analysis depend on an assumption of the applicability of the spectral model in Eq. (5), while the two-frequency results are model independent.

The multi-frequency analysis shows no evidence for a loss of correlation, within the limits provided by the analysis of the E2E simulations. As discussed, these new limits are much tighter than those obtained from the 217- and 353 GHz correlation ratio in Appendix B. However, current limits are still consistent with (i.e. still allow the presence of) significant variations of the dust spectral index over the sky. To illustrate this statement quantitatively, we have computed for the noise-free FFP10 dust polarization maps (Appendix A.1), built from 353 GHz polarization templates computed using the Vansyngel et al. (2017) model. These 353 GHz templates were scaled to other frequencies using maps of dust temperature and spectral index that were derived from the analysis of dust total intensity maps in Planck Collaboration Int. XLVIII (2016) and Planck Collaboration XI (2014). The standard deviations of the dust spectral index for our six sky regions, measured using the 217–353 GHz colour ratio of model maps smoothed to a 1° resolution, are in the range σ(β_d) = 0.092 ± 0.005. Nevertheless, the values of that we obtained, listed in Table 5, are within the lower limits inferred from the E2E simulations.

Frequency decorrelation might result from variations of the spectral index both across the sky and along the line of sight. In the FFP10 maps, only the former is taken into account, and thus polarization angles do not vary with frequency (Tassis & Pavlidou 2015, PL). We also note again that the results of our multi-frequency analysis depend in detail on the adopted spectral model (Eq. (4)).

We have also used the 300 E2E realizations to compute the cross-correlation of measured for our six sky regions (Fig. 19). As found for the correlation analysis between the 217- and 353 GHz data discussed in Appendix B and by Sheehy & Slosar (2018), the results from the multi-frequency fit are also correlated between sky regions, which makes sense, of course, because they are nested.

Fig. 19. Cross-correlation factor between sky regions, determined from the multi-frequency fits expressed as percentages. As found for the correlation analysis between the 217- and 353 GHz data presented in Appendix B, the results from the multi-frequency fit are correlated between the six nested sky regions.

6.2. Perspective for on-going and future CMB experiments

Here, we discuss the implications of our multi-frequency analysis of Planck dust polarization for on-going and future CMB experiments that are designed to search for primordial B modes. A somewhat comforting view concerning the complexity of dust polarization as a CMB foreground is suggested by two of our results. First, the data show no departure from a one-parameter MBB emission law (with a single fixed temperature) for the dust polarization SED spanning from 353 GHz to below 70 GHz. Second, the data do not provide evidence for frequency decorrelation.

For our largest sky region, LR71, our lower limit on the correlation ratio is tightest and quite close to unity for the multipole range relevant to the search for primordial B modes at the recombination peak. Using Eq. (5), this limit translates to limits on the correlation ratio between the frequencies (95, 150, 220 GHz) used in ground-based and balloon-borne experiments, namely 0.996 and 0.979 for and , respectively.

Sub-orbital CMB experiments typically target smaller sky regions, for which frequency decorrelation is less constrained by the Planck data. If the LR71 Planck limit applies to these cleaner sky regions, then frequency decorrelation of dust polarization might not be a problem for CMB experiments aiming at a primordial B-mode detection limit of r ≃ 0.01 at the recombination peak. To quantify this, we consider LR24 as representative of clean sky regions used by sub-orbital experiments for their CMB studies. For this region, at ℓ = 80, and for 150 and 95 GHz, the dust power is 45 and 7 times the primordial B-mode signal for r = 0.01, respectively. Combining these factors with the corresponding values of and for LR71, the potential bias from decorrelation when combining data at 95, 150, and 220 GHz could be smaller than the CMB B-mode signal for r = 0.01.

We caution that this view might be too optimistic because the frequency decorrelation might not be homogeneous over the sky. Indeed, the tight limit that Planck data provide for LR71 might not apply to smaller sky regions. In particular, decorrelation could have a large statistical variance for small sky regions, where there is not a large amount of averaging over relevant scales in the ISM. Furthermore, the Planck limits in Table 5 are dependent on the adopted model and so frequency decorrelation of dust polarization at microwave frequencies could be larger than we have estimated.

It is worth stressing that the Planck sensitivity precludes identifying how difficult the component separation will be for more ambitious experiments – CMB Stage IV (Abazajian et al. 2016) and a future space mission, for example LiteBIRD (Ishino et al. 2016) – that are targeting a B-mode detection limit r ≃ 0.001 using the recombination and reionization peaks. Given the limitations of available data, it is essential to continue to make progress in assessing the component-separation problem by using increasingly realistic models developed in relation to the astrophysics of the dusty magnetized ISM.

7. Conclusions

Using the Planck PR3 polarization maps, this paper has extended the characterization of Galactic dust polarized emission that is foreground to CMB polarization. Our data analysis is validated using E2E simulations, where the mapmaking pipeline is run on simulated data derived by combining fixed maps of polarized sky emission with independent realizations of the data noise and systematics. This final section summarizes the main results of our study.

The power spectra of dust polarized emission (EE, BB, TE, TB, and EB) are measured for six nested high-Galactic latitude regions, covering a range of sky fractions from 24 to 71%, over the multipole range (ℓ ≤ 600) relevant to the analysis of E- and B-mode CMB polarization associated with reionization, recombination, and lensing. We present power-law fits to the angular power spectra that reveal statistically significant variations of the exponents over sky regions and a difference between the values for the EE and BB spectra, which for the largest sky region are α_EE = −2.42 ± 0.02 and α_BB = −2.54 ± 0.02, respectively. The difference persists in the weighted mean values, α_EE = −2.39 and α_BB = −2.51. The small difference between the two exponents is not unexpected because the filamentary structures in the cold neutral interstellar medium have mainly E-mode polarization, due to the statistical alignment of the magnetic field orientation with matter. The BB power scales as the square of the mean intensity of the region, ⟨I₃₅₃⟩. The spectra show that the TE correlation and the E/B power asymmetry discovered by Planck extend to low multipoles, which were not included in the earlier Planck polarization papers due to residual data systematics. The weighted mean value of TE/EE is 2.76 ± 0.05. The mean TE correlation ratio is  = 0.357 ± 0.003, with a scatter of around 0.1 between measured values, but no systematic dependence on multipole down to the lowest ℓ bins or on sky region. We also report a significant TB signal with a TB/TE ratio of approximately 0.1 and correlation ratio about 0.05.

Combining data from Planck and WMAP, we characterize the mean SED of polarized Galactic foregrounds for the six sky regions as a function of multipole, for ℓ <  160. Our spectral model takes into account polarized synchrotron emission and its correlation with polarized dust emission. The results of this analysis quantify the challenge of the component-separation procedure required for measuring the low-ℓ CMB E-mode reionization signal and detecting the reionization and recombination peaks of primordial CMB B modes.

In our analysis, we do not find systematic variations of the polarized dust SED with multipole value or with sky region. The mean dust spectral index is  = 1.53 ± 0.02 for a dust temperature of 19.6 K. The systematic error follows from the uncertainties on the polarization efficiencies of HFI and includes the uncertainty from the CMB subtraction. The dust SED in polarization from blind component separation is remarkably well fit by a single temperature MBB emission law from 353 to 44 GHz with a similar index. The difference between the indices for polarization and total intensity is small and not of high statistical significance, . This result suggests that the emission from a single grain type dominates the long-wavelength emission in both polarization and total intensity. It constrains dust models involving multiple dust components (e.g. separate carbon and silicate grains), magnetic dipole emission and variations of the spectral index of the dust emissivity with frequency. Detailed modelling, beyond the scope of this paper, is required to quantify these constraints.

We analyse the correlation of the dust polarization maps across microwave frequencies by fitting cross-spectra between Planck data at 100, 143, 217, and 353 GHz. We find no evidence for a loss of correlation with frequency, within limits provided by the analysis of the E2E simulations. These new results provide tighter lower limits to the correlation ratio than we obtain from the comparison of the 217- and 353 GHz maps alone. If the Planck limit on decorrelation for the largest sky region applies to the smaller sky regions observed by sub-orbital experiments, then frequency decorrelation of dust polarization might not be a problem for CMB experiments aiming at a primordial B-mode detection limit at the r ≃ 0.01 level using the recombination peak. However, the sensitivity of Planck prevents us drawing any conclusions about how difficult the component-separation problem will be for more ambitious experiments targeting lower levels for r.

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration. This research has received funding by the Agence Nationale de la Recherche (ANR-17-CE31-0022).

References

Appendix A: Data simulations and uncertainties

In this appendix, we detail how the uncertainties in the data are assessed and propagated to the power spectra used in Sects. 3 and 4. End-to-end (E2E) simulations are used throughout the paper in order to check for potential biases in our data analysis procedure and evaluate uncertainties in our results.

A.1. Maps of data uncertainties

We use E2E simulations to build uncertainty maps for the PR3 HFI polarization data, which include both residual systematics and instrumental noise. Simulations of Planck timelines are computed by combining the scanning strategy with models of known systematic effects and models of the sky emission (Appendix A.2). The systematic effects included in the simulations are described in Planck Collaboration III (2020). One single sky model, referred to as the FFP10 sky model in this paper, is used. For all of the sky components except dust polarization, in particular polarized synchrotron emission, we used the latest version of the Planck sky model described in Planck Collaboration XII (2016).

The dust model maps are built as follows. The Stokes I map at 353 GHz is the dust total intensity Planck map obtained by applying the Generalized Needlet Internal Linear Combination (GNILC) method of Remazeilles et al. (2011) to the 2015 release of Planck HFI maps (PR2), as described in Planck Collaboration Int. XLVIII (2016), and subtracting the monopole of the cosmic infrared background (Planck Collaboration VIII 2016). For the Stokes Q and U maps at 353 GHz, we started with one realization of the statistical model of Vansyngel et al. (2017). The portions of the simulated Stokes Q and U maps near Galactic plane were replaced by the Planck 353 GHz PR3 data. The transition between data and the simulations was made using a Galactic mask with a 5° apodisation⁷, which leaves 68% of the sky unmasked at high latitude. Furthermore, on the full sky, the large angular scales in the simulated Stokes Q and U maps were replaced by the Planck data. Specifically, the first ten multipoles were the Planck data, while over ℓ = 10 to ℓ = 20 the simulations were introduced smoothly using the function (1 + sin [π(15−ℓ)/10])/2. The resulting 353 GHz Stokes maps were scaled to other frequencies using the maps of dust temperature and spectral index, coming from fitting the SED of dust total intensity in Planck Collaboration Int. XLVIII (2016) and Planck Collaboration XI (2014), respectively, which introduces significant spectral variations as discussed in Sect. 6.1. Hereafter, we refer to these simulated maps as the FFP10 dust polarization maps.

Independent realizations of the detector noise and data systematics are computed for each simulation of the HFI timelines, keeping the sky emission components (including the CMB) the same. The mapmaking pipeline is run on the simulated timelines to produce simulated maps. This process is repeated to obtain 300 realizations of simulated maps at each of the four HFI polarized channels, at 100, 143, 217, and 353 GHz. We compute difference maps by subtracting the input sky model from each of the simulated maps. The 300 difference maps are independent E2E realizations of HFI uncertainty maps, including both detector noise and systematic effects.

For illustration, one difference map obtained for one random realization of the E2E simulations at 353 GHz is presented in Fig. A.1. Figure A.2 presents the mean EE and BB power spectra of the residual maps, computed over the set of 300 simulations. These spectra are compared with those computed with the half-difference of the two half-mission PR3 maps, which quantify the instrumental noise in the data. The noise spectra computed from the half-differences between half-mission and odd-even survey maps are very similar, as shown in Fig. A.3. The excess power in the E2E spectra at low multipoles corresponds to residual uncorrected systematics, which are taken into account in the HFI E2E uncertainty maps.

Fig. A.1. Difference between the Stokes Q and U output maps and the sky model inputs at 1° resolution, for one E2E simulation at 353 GHz. Such pairs of difference maps are used to quantify statistical and systematic uncertainties in our analysis of Planck HFI data.

Fig. A.2. EE (left panels) and BB (right panels) power spectra from the HFI uncertainty maps described in Appendix A.1. The data points and error bars represent the mean and standard deviation of the power in each multipole bin computed over the 300 realizations for the LR42 (red) and LR71 (blue) sky regions. The dashed lines represent analytical fits (a power law plus a constant) to the data points. For comparison, the plots also present the spectra of the difference between half-mission maps (solid lines, labelled “HMD”).

Fig. A.3. Like Fig. A.2, but comparing EE and BB power spectra of the data noise estimated from the half-differences between odd and even survey maps (OED data points and dashed lines) with that estimated from the half-mission maps (solid lines and HMD data points, as in Fig. A.2).

We use LFI and WMAP maps to separate dust and synchrotron polarized emission and quantify the correlation between the two sources of emission. Because E2E simulations are not available for these data, we compute maps of uncertainties from Gaussian realizations of the data noise. Power spectra of this noise are derived from the half-difference of half-mission Planck LFI maps and the difference of year maps for WMAP. We note that it is easy to produce a large number (1000 or more) of data realizations with Gaussian data noise, while only 300 E2E realizations are available for HFI.

A.2. E2E simulated maps

The uncertainty maps for all four HFI frequencies described in the previous section are combined with a simple model of dust and synchrotron polarized emission to produce what we call E2E simulated maps. We use these to quantify error bars on the power spectra and validate our data analysis. The same single set of 300 E2E simulated maps is used for all sky regions throughout the paper.

A single sky model is used for all simulations. Stokes maps of the dust polarization at 353 GHz are computed using the model of Vansyngel et al. (2017). These model maps are scaled to the other frequencies with a single SED for all sky pixels, namely an MBB emission law with a spectral index of 1.59 and a temperature of 19.6 K based on the Planck data analysis in Planck Collaboration Int. XVII (2014), Planck Collaboration Int. XXII (2015). This model, unlike the FFP10 dust polarization maps (Appendix A.1), has no spectral variations. The spatial template for the synchrotron polarized emission is derived from the Planck sky model, as in Planck Collaboration XII (2016). The SED of the synchrotron polarized emission is a power law with a spectral index of −3.0 for all sky pixels. Independent realizations of the CMB polarization maps, computed from the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016), are added to the Stokes maps of the dust and synchrotron emission. By combining the uncertainty maps from the E2E simulations (Appendix A.1) with Galactic polarization maps and CMB realizations distinct from those used in the simulations (we note that the Stokes I dust map is unchanged), we erase potential correlations between data systematics and the polarized sky emission. Such correlations have been shown to have negligible impact on the CMB data analysis (Planck Collaboration Int. XLVI 2016; Planck Collaboration III 2020). We also checked that the correlation between dust polarization and residual systematics is not a dominant uncertainty at 353 GHz.

A.3. Uncertainties propagated to power spectra in the data analysis

We compute EE and BB power spectra of a number of E2E simulated maps, which combine the sky model (including CMB polarization) with independent realizations of the statistical Gaussian noise for LFI and WMAP or the E2E uncertainty maps for HFI. The spectra are computed with XPol (Tristram et al. 2005) for the same multipole bins and sky regions used for the data analysis. The CMB contribution is subtracted from the power spectra using the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016).

The dispersion of values computed over a set of 300 realizations define the statistical uncertainties that we adopt for the corresponding 𝒞_ℓ coefficient measured on the data maps. Because each set of Stokes I, Q, and U E2E simulated maps includes a different realization of the CMB polarization, the uncertainties include the cosmic variance of the CMB. These uncertainties are used in the power-law fits in Sect. 3.2 and the fits to the spectral model in Sect. 4.2.

The data simulations are also used to check for potential biases in the fit of the spectral model. For example, the results presented in Fig. A.4 for the parameters of the spectral fit to the E2E simulated HFI maps show that for all multipole bins we recover the input spectral indices of dust and synchrotron polarized emission (β_d and β_s) of the sky model without any bias. We point out that this validation has been done using the least-squares MPFIT fitting routine and not the more computationally-intensive MCMC code used in Sect. 4.2 to fit the data. We have checked that the two methods produce consistent determinations of the model parameters when applied to the same data.

Fig. A.4. Results of the spectral fit to simulated data with uncertainties derived from the E2E simulated maps. Distinct symbols and colours represent the different sky regions (see top row). The dashed lines in plots of βs and βd are the input values of our sky model. We retrieve these input values without any bias. The sky model has a non-zero correlation between dust and synchrotron polarization at low ℓ, which is consistent with the value of ρ found by fitting the simulated data.

Appendix B: Two-frequency analysis of the spectral correlation ratio between 217 and 353 GHz

PL computed the spectral correlation ratio , defined as

(B.1)

If the B-mode emission is perfectly correlated between the two frequencies, then , whereas a value lower than 1 is indicative of a correlation that is only partial. PL interpreted their results as evidence for decorrelation and spatial variations of dust polarization between 217 and 353 GHz, over multipoles relevant to the search for primordial B modes at the recombination peak. In the course of our analysis of the PR3 data we found that this conclusion was compromized as described below, so that the significance of the decorrelation was overstated. Sheehy & Slosar (2018) discovered this independently.

In this Appendix, using the new Planck maps, we update and extend the PL analysis (Appendix B.1). The E2E simulations are used to assess the uncertainty of and the statistical significance of the results (Appendix B.2). We note that does not depend on the 1.5% uncertainty on the 353 GHz polarization efficiency.

B.1. Measured ratios on PR3 polarizations maps

Here, we compute using the PR3 Planck data for five broad ranges of multipoles, namely ℓ = 4–11, 11–50, 50–160, 160–320, and 320–500. The three last ℓ bins are common to the analysis reported by PL using the PR2 Planck data. The lowest two ℓ bins are new.

The values of are listed for the six sky regions, LR24 to LR71, in Table C.5. In Fig. B.1, the ratios measured for the first four ℓ bins are plotted versus the mean hydrogen column density, N_H. The bottom left panel, for the bin ℓ = 50–160, is directly comparable to the corresponding plot from the PL analysis of the PR2 Planck data in their Fig. 3. For this common ℓ bin, we find results consistent with the earlier analysis in PL, that the departure from 1 (i.e. the apparent spectral decorrelation) increases with decreasing column density. As in PL, we also find that the apparent spectral decorrelation increases with increasing multipole, as illustrated in Fig. B.2. However, the two new lowest ℓ bins show very little decorrelation.

Fig. B.1. Spectral correlation ratio versus mean hydrogen column density for the six sky regions. Each panel presents the results for one of the ℓ bins, ℓ = 4–11, 11–50, 50–160, and 160–320. The dark and light grey bars represent the 68% and 95% confidence intervals computed over the 300 E2E simulations, while the blue segments mark the median values. The bottom left panel for the ℓ = 50–160 bin is directly comparable to the corresponding plot from the PL analysis, their Fig. 3.

B.2. Statistical significance of the Planck results

We have used the E2E simulated HFI maps described in Appendix A.2 to compute probability distribution functions of the ratio and thereby quantify the statistical significance of the values found from the Planck data. The histograms are shown for the five bins, ℓ = 4–11, 11–50, 50–160, 160–320, and 320–500, in Fig. B.3. The dust polarization model used in the simulations has the same SED for each sky pixel, i.e. a perfect correlation. The plots in Fig. B.1 show the 68% and 95% confidence intervals for the values measured on the 300 E2E simulations. The shift of the median values, marked by the dashed lines from , is a bias due to data uncertainties. The bias and the width of the distribution are both larger in this analysis (based on the E2E simulations) than in PL (based on uncertainty maps derived from Gaussian noise realizations). These differences result from the fact that E2E simulations take into account uncorrected systematics, ignored in Gaussian noise realizations. For our two lowest ℓ bins, the data values are within the 68% interval of the simulation results. For the ℓ = 50–160 and 160–320 bins, the values measured on the Planck data are intermediate between the 68% and 98% intervals of the simulation results, i.e. between the 1 and 2 σ limits for a Gaussian distribution⁸.

Fig. B.2. Spectral correlation ratio versus multipole for the six sky regions. Like in PL, we find that the apparent decorrelation of the dust polarization measured by the difference increases for decreasing , that is towards low-brightness regions at high Galactic latitude.

The E2E simulations show that for a given ℓ bin the values are highly correlated over sky regions. The measured cross-correlation ratios between regions are displayed in Fig. B.4 for all our ℓ bins. Values range from 50% to more than 90%. This correlation follows from the fact that the sky regions we use are nested. The data points in Fig. B.1 are not independent, a fact also identified by Sheehy & Slosar (2018).

Fig. B.3. Distribution of the correlation ratios on the six sky regions for each of the five ℓ bins. The histograms are computed from the 300 E2E simulations using half-mission data splits. Histograms computed using odd-even surveys give similar results. The values derived from the data are displayed as vertical solid lines. The dashed lines represent the median value for the simulations.

The values in independent ℓ bins are uncorrelated, but noise introduces a systematic trend, where decreases for increasing multipole, which accounts for the systematic trend in Fig. B.2.

To illustrate the correlation between nested regions, we have measured on two independent sets of sky regions, corresponding to the northern and southern Galactic parts of the LR52, LR62, and LR71 regions. The results of this analysis are presented in Fig. B.5. The cross-correlation factors between these north-south splits, computed from the E2E simulations, are plotted in Fig. B.6. The data points within a given set (north or south) are correlated, but the two sets are independent. For example, for the ℓ = 50–160 bin, the values of for the northern regions are consistently lower than those for the corresponding southern regions.

Fig. B.4. cross-correlation factor between sky regions for three ℓ bins, expressed as percentages.

Fig. B.5. Spectral correlation ratio versus the mean hydrogen column density for northern and southern splits of three large sky regions. The dark and light grey bars represent the 68% and 95% intervals computed over the 300 E2E simulations, while the blue segments mark the median values.

Fig. B.6. cross-correlation measured for the north-south splits of the L52, L62, and L71 sky regions for three ℓ bins, expressed as percentages.

In Table C.5, for each ℓ bin and sky region we list the probability, labeled PTS, to obtain correlation ratios smaller than the Planck measurements, based on the 300 E2E simulations. The lower the value of the PTS, the more significant is the measurement. Table C.6 lists the same probability for the north-south splits.

As in Sheehy & Slosar (2018), we conclude that our Planck spectral correlation analysis of the 217- and 353 GHz polarization maps does not provide statistically significant evidence for a departure of the spectral correlation ratios from a perfect correlation (i.e. a ratio of 1). Lower limits are listed for our set of sky regions and ℓ bins in Table B.1. These limits are from the 2.5th percentile of the 300 E2E simulations (relating to the 95% confidence interval), i.e. the values of are smaller than these limits for only 7.5 realizations out of the 300. The values of for the ℓ bin 50–160 presented in Fig. 3 of PL are within these limits.

Appendix C: Additional tables

All Tables

All Figures

	Fig. 1. All-sky map showing the sky regions used to measure power spectra, indicated with colours varying from yellow to orange and dark-red. The white region represents the area where the CO line brightness is larger than 0.4 K km s−1, which is excluded from all the sky regions in our analysis. The blue dots represent the areas masked around point sources.
In the text

	Fig. 2. CMB-corrected EE (red diamonds), BB (blue squares), and TE (black circles) power spectra at 353 GHz, for each of the six sky regions that we analyse. The dashed lines represent power-law fits to the data points from ℓ = 40 to 600. The exponents of these fits, αTE, αEE, and αBB, appear on each panel.
In the text

	Fig. 3. Power spectra, as in Fig. 2, but for the northern and southern parts of the LR42, LR52, LR62, and LR71 regions.
In the text

Fig. 4. Scaling of the BB power at ℓ = 80 versus the mean dust total intensity at 353 GHz. The dashed black line is a power-law fit to values for the six sky regions in our analysis (blue dots) with an exponent of two. Also shown are the values for the N–S splits of the regions in Fig. 3 (blue triangles). These results are complemented by the measurement (red diamond) over the southern Galactic cap (fsky = 8.5%) by Ghosh et al. (2017) and that for the BICEP field (black square) after BICEP2 and Keck Array Collaborations (2016).

In the text

	Fig. 5. TE correlation ratio versus multipole. The data points are plotted using distinct symbols and colours (see legend at the top) for each of the six sky regions. The error bars are derived from the E2E simulations.
In the text

	Fig. 6. Power spectra of TB (red diamonds) and EB (blue squares) at 353 GHz for the six sky regions. The error bars are derived from the E2E simulations. A power-law fit to the TB data (solid red line) reveals an overall positive TB signal, not seen in the E2E simulations. The EB power (solid blue line fit) is consistent with zero.
In the text

Fig. 7. BB cross-spectra 𝒟ℓ(ν1 × ν2) versus the effective frequency νeff = (ν1 × ν2)0.5, for the LR62 sky region and two multipole bins: ℓ = 4-11 (top plot) and 40–59 (bottom). Yellow and blue colours represent data values from single and inter-frequency cross spectra, respectively. Bottom panel: within each plot shows the residuals from the fits normalized to the 1 σ uncertainty of each data point. Lower frequency data (left) points are dominated by the SED of synchrotron polarized emission, while higher frequency (right) data characterize dust polarized emission, and those at the centre characterize the correlation between the two sources of emission. Differences between the two plots illustrate that both the ratio between synchrotron and dust power and the correlation between these two sources of polarized emission decrease for increasing multipoles.

In the text

Fig. 8. Posterior distribution for each of the parameters of the spectral model in Eq. (2), as obtained through the MCMC fitting algorithm for BB data points. The MCMC results illustrated here are for the LR62 region and the multipole bin ℓ = 40–59, one of the two cases shown in Fig. 7. The diagonal shows the probability distribution of each parameter. Median values are As = 0.6 ± 0.1, Ad = 137 ± 2, βs = −3.15 ± 0.17, βd = 1.50 ± 0.02, and ρ = 0.17 ± 0.04.

In the text

Fig. 9. Amplitudes of EE and BB power spectra for dust and synchrotron emission at 353 and 30 GHz, respectively, shown for each sky region and each multipole bin. The As and Ad parameters of our spectral model from Eq. (2) are thermodynamic (CMB) temperature in μK2. Where the synchrotron amplitude is compatible with zero at the 1 σ level, we report an upper limit on As (68% confidence limit) with triangles pointing down.

In the text

	Fig. 10. Fit parameters ρ and βd for E- and B-mode polarization versus multipole. Open symbols for ρ represent the cases where the synchrotron amplitude is compatible with zero, making it difficult to measure the correlation.
In the text

Fig. 11. Dust and synchrotron E-mode power versus multipole. The dust power at 95 and 150 GHz and that of synchrotron at 95 GHz are compared with the CMB E-mode signal (red-line) computed for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016) and a Thompson scattering optical depth τ = 0.055 from Planck Collaboration Int. XLVI (2016). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018).

In the text

Fig. 12. Dust and synchrotron B-mode power versus multipole. The dust power at 95 and 150 GHz, and that of synchrotron at 95 GHz are compared with CMB B modes from primordial gravitational waves (grey lines) for three values of the tensor-to-scalar ratio, r = 0.1, 0.01, and 0.001, and from lensing (blue line) for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018).

In the text

Fig. 13. Dust and synchrotron B-mode power versus frequency for two multipole bins: ℓ = 4–11 (top) and 60–79 (bottom). The coloured bands show the range of power measured from the smallest (LR24) to the largest (LR71) sky regions in our analysis. The lower limit of the synchrotron band is derived from the S-PASS data analysis in Krachmalnicoff et al. (2018). The primordial CMB B-mode signal, averaged within the appropriate ℓ bin, is plotted with dashed lines for three values of the tensor-to-scalar ratio: r = 0.1; 10−2; and 10−3. The solid line represents the lensing B-mode signal for the Planck 2015 ΛCDM model (Planck Collaboration XIII 2016).

In the text

Fig. 14. Dust SEDs for E- and B-mode polarization derived from the SMICA component-separation procedure (Planck Collaboration IV 2020). The two grey lines represent MBB fits to the E- (red diamonds) and B-mode (blue squares) data points with a temperature of 19.6 K. The polarization spectral index derived from the fits is  = 1.53 ± 0.02. The residuals to each fit, normalized to the 1 σ data uncertainty, are plotted in the lower panel.

In the text

Fig. 15. Comparison of spectral indices of dust polarized emission and total intensity. The spectral indices are derived from the 353-to-217 GHz colour ratio. Plots to the left show the results obtained from our simulated maps, and the ones to the right are from the Planck data. Distinct symbols are used to represent each of the six sky regions, as in Fig. 10. For the simulations, the dashed lines represent the input dust spectral indices (, ). For the data, the dashed lines represent the mean measured dust spectral indices (, ).

In the text

Fig. 16. Illustration of the consistency between our analysis and component separation methods. The probability distribution of the 217–353 GHz colour ratio for dust polarized intensity, computed over the LR71 sky region from Planck CMB-subtracted maps smoothed to a 3° beam, is plotted for each of the component separation methods in Planck Collaboration IV (2020). The vertical line is the value derived from our analysis. For the unit conversion factors and color corrections, and our modification of the 217 GHz polarization efficiency, it corresponds to the spectral index from Sects. 5.1 and 5.2. The width of the line represents the error bar.

In the text

	Fig. 17. Posterior distribution for each of the parameters of the spectral model with decorrelation given in Eq. (5), as obtained through the MCMC fitting algorithm for BB data points. The MCMC results illustrated here are for the LR62 region and the multipole range 50–160. Median values are Ad = 97.1 ± 1.2, βd = 1.54 ± 0.02, and .
In the text

Fig. 18. Distribution of the correlation ratios inferred from δd on the six sky regions for the ℓ range 50–160. The histograms are computed from the 300 E2E simulations using half-mission data splits. The dashed lines represent the median values on each sky region. This median value, μ, and the standard deviation, σ, are printed in the upper right of each panel. The lower limits on in Table 5 are derived from the 2.5th percentile of the distribution for each sky region.

In the text

	Fig. 19. Cross-correlation factor between sky regions, determined from the multi-frequency fits expressed as percentages. As found for the correlation analysis between the 217- and 353 GHz data presented in Appendix B, the results from the multi-frequency fit are correlated between the six nested sky regions.
In the text

	Fig. A.1. Difference between the Stokes Q and U output maps and the sky model inputs at 1° resolution, for one E2E simulation at 353 GHz. Such pairs of difference maps are used to quantify statistical and systematic uncertainties in our analysis of Planck HFI data.
In the text

Fig. A.2. EE (left panels) and BB (right panels) power spectra from the HFI uncertainty maps described in Appendix A.1. The data points and error bars represent the mean and standard deviation of the power in each multipole bin computed over the 300 realizations for the LR42 (red) and LR71 (blue) sky regions. The dashed lines represent analytical fits (a power law plus a constant) to the data points. For comparison, the plots also present the spectra of the difference between half-mission maps (solid lines, labelled “HMD”).

In the text

	Fig. A.3. Like Fig. A.2, but comparing EE and BB power spectra of the data noise estimated from the half-differences between odd and even survey maps (OED data points and dashed lines) with that estimated from the half-mission maps (solid lines and HMD data points, as in Fig. A.2).
In the text

Fig. A.4. Results of the spectral fit to simulated data with uncertainties derived from the E2E simulated maps. Distinct symbols and colours represent the different sky regions (see top row). The dashed lines in plots of βs and βd are the input values of our sky model. We retrieve these input values without any bias. The sky model has a non-zero correlation between dust and synchrotron polarization at low ℓ, which is consistent with the value of ρ found by fitting the simulated data.

In the text

Fig. B.1. Spectral correlation ratio versus mean hydrogen column density for the six sky regions. Each panel presents the results for one of the ℓ bins, ℓ = 4–11, 11–50, 50–160, and 160–320. The dark and light grey bars represent the 68% and 95% confidence intervals computed over the 300 E2E simulations, while the blue segments mark the median values. The bottom left panel for the ℓ = 50–160 bin is directly comparable to the corresponding plot from the PL analysis, their Fig. 3.

In the text

	Fig. B.2. Spectral correlation ratio versus multipole for the six sky regions. Like in PL, we find that the apparent decorrelation of the dust polarization measured by the difference increases for decreasing , that is towards low-brightness regions at high Galactic latitude.
In the text

	Fig. B.3. Distribution of the correlation ratios on the six sky regions for each of the five ℓ bins. The histograms are computed from the 300 E2E simulations using half-mission data splits. Histograms computed using odd-even surveys give similar results. The values derived from the data are displayed as vertical solid lines. The dashed lines represent the median value for the simulations.
In the text

	Fig. B.4. cross-correlation factor between sky regions for three ℓ bins, expressed as percentages.
In the text

	Fig. B.5. Spectral correlation ratio versus the mean hydrogen column density for northern and southern splits of three large sky regions. The dark and light grey bars represent the 68% and 95% intervals computed over the 300 E2E simulations, while the blue segments mark the median values.
In the text

	Fig. B.6. cross-correlation measured for the north-south splits of the L52, L62, and L71 sky regions for three ℓ bins, expressed as percentages.
In the text