# Dust polarization spectral dependence from Planck HFI data. Turning point on CMB polarization foregrounds modelling

Alessia Ritacco, François Boulanger, Vincent Guillet, Jean-Marc Delouis, Jean-Loup Puget, Jonathan Aumont, Léo Vacher

► To cite this version:

Alessia Ritacco, François Boulanger, Vincent Guillet, Jean-Marc Delouis, Jean-Loup Puget, et al. Dust polarization spectral dependence from Planck HFI data. Turning point on CMB polarization foregrounds modelling. Astronomy and Astrophysics - A&A, 2023, 670, pp.A163. 10.1051/0004-6361/202244269. hal-03710233

HAL Id: hal-03710233

https://hal.science/hal-03710233

Submitted on 10 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dust polarization spectral dependence from Planck HFI data

## Turning point for cosmic microwave background polarization-foreground modeling

Alessia Ritacco<sup>1,2,3</sup>①, François Boulanger<sup>2</sup>①, Vincent Guillet<sup>3,4</sup>①, Jean-Marc Delouis<sup>5</sup>①, Jean-Loup Puget<sup>2,3</sup>, Jonathan Aumont<sup>6</sup>①, and Léo Vacher<sup>6</sup>①

1 INAF-Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius, Italy
e-mail: alessia.ritacco@inaf.it

2 Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

3 Institut d'Astrophysique Spatiale, CNRS, Université Paris-Saclay, CNRS, Bât. 121, 91405 Orsay, France

4 Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

5 Laboratoire d'Océanographie Physique et Spatiale (LOPS), Univ. Brest, CNRS, Ifremer, IRD, Technopole Brest-Iroise, 29280 Plouzane, France

6 IRAP, Université de Toulouse, CNRS, CNES, UPS, 9 Av. du Colonel Roche, 31400 Toulouse, France

Received 14 June 2022 / Accepted 11 December 2022

## ABSTRACT

The search for the primordial *B*-modes of the cosmic microwave background (CMB) relies on the separation from the brighter foreground dust signal. In this context, the characterization of the spectral energy distribution (SED) of thermal dust in polarization has become a critical subject of study. We present a power-spectra analysis of *Planck* data, which improves upon previous studies by using the newly released SRo112 maps that include corrections on residual data systematics and by extending the analysis to regions near the Galactic plane. Our analysis focuses on the lowest multipoles between $l = 4$ and 32, as well as three sky areas with sky fractions of $f_{sky} = 80\%, 90\%,$ and $97\%$. The mean dust SED for polarization and the 353 GHz *Q* and *U* maps are used to compute residual maps at 100, 143, and 217 GHz, highlighting variations of the dust polarization SED on the sky and along the line of sight. Residuals are detected at the three frequencies for the three sky areas. We show that models based on total-intensity data end up underestimating (by a significant factor) the complexity of dust polarized CMB foreground. Our analysis emphasizes the need to include variations of the polarization angles of the dust polarized CMB foreground. The frequency dependence of the *EE* and *BB* power spectra of the residual maps yields further insight. We find that the moments expansion to the first order of the modified black-body (MBB) spectrum provides a good fit to the *EE* power-spectra. This result suggests that the residuals could follow mainly from variations of the dust MBB spectral parameters. However, this conclusion is challenged by cross-spectra showing that the residuals maps at the three frequencies are not fully correlated, as well as the fact that the *BB* power-spectra do not match the first order moment expansion of a MBB SED. This work sets new requirements for simulations of the dust-polarized foreground and component separation methods, showing that a significant refinement to the dust modeling is necessary to ensure an unbiased detection of the CMB primordial *B*-modes at the precision required by future CMB experiments. Further works would also be required to theoretically model the impact of polarization-angle variations on the *EE* and *BB* power spectra of residual maps.

**Key words.** cosmic background radiation – cosmology: observations – diffuse radiation – polarization

## 1. Introduction

One of the outstanding questions in cosmology concerns the existence of primordial gravitational waves, as predicted by the theory of cosmic inflation (Guth 1981; Linde 1982). Although the primordial gravitational waves are not directly detectable with foreseen experiments, they are expected to leave an imprint that is seen as a curl-like pattern in the cosmic microwave background (CMB) polarization anisotropies, referred to as primordial *B*-modes, which can be measured. This is the main goal of present and future CMB experiments, including the LiteBIRD satellite (LiteBIRD Collaboration 2022) and BICEP/Keck (Ade et al. 2022), the Simons Observatory (Ade et al. 2019), and CMB-Stage 4 (Abazajian et al. 2016) from the ground. These experiments are aimed at measuring the tensor-to-scalar ratio, *r*, which represents the amplitude of the primordial tensor (*B*-modes) relative to scalar perturbations (*E*-modes)

of the CMB. This parameter, related to the energy scale of inflation, is expected to be in the range between $10^{-2}$ to $10^{-4}$ (Kamionkowski & Kovetz 2016; LiteBIRD Collaboration 2022). A very high control and subtraction of instrumental systematic effects and Galactic foregrounds is required for an unbiased measurement of such low *r* values.

The dust emission represents a major obstacle because its amplitude is much higher than that of the primordial *B*-modes signal (Planck Collaboration XI 2020). In this context, the characterization of the spectral energy distribution (SED) of dust has become a critical step in the search for *B*-modes. The *Planck* data have shown that the mean dust SED for polarization and total intensity are very close and that they are both well fitted by a modified black body (MBB) law (Planck Collaboration Int. XXII 2015; Planck Collaboration XI 2020). In addition, far-IR polarization measurements obtained by the BLASTPol balloon-borne experiment (e.g., Ashton et al. 2018), have shown that the
dust polarization fraction is roughly constant between 250 µm and 3 mm (100 GHz). These remarkable results suggest that the emission from a single grain type dominates the long-wavelength emission in both polarization and total intensity (Guillet et al. 2018; Hensley & Draine 2022).

Maps of MBB parameters (dust spectral index and temperature) have been obtained by fitting the total intensity Planck data (e.g., Planck Collaboration X 2016; Planck Collaboration Int. XLVIII 2016). The Planck 2018 data release (hereafter, PR3, Planck Collaboration III 2020) does not provide comparable constraints on MBB parameters in terms of the polarization due to an insufficient signal-to-noise ratio (S/N) and instrumental systematics (Osumi et al. 2021). These limiting factors are emphasized by the lack of polarization maps at 545 and 857 GHz which would be useful in constraining the dust SED. The frequency dependence of dust polarization data also involves polarization angles. When the dust SED and magnetic field orientation vary within the beam, the frequency scaling of the Stokes Q and U parameters may differ (Ichiki et al. 2019; Vacher et al. 2023). Integration along the line of sight may therefore induce variations in the polarization angle with frequency (Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017). Pelgrims et al. (2021) provided first observational evidence of this effect analyzing Planck data toward the lines of sight with multiple velocity components in H I emission. Additional emission components (e.g., magnetic dipole and CO emission Puglisi et al. 2017) can also contribute to SED variations (Hensley & Bull 2018).

Any SED variations on the sky and along the line of sight induce a decorrelation between dust emission at different frequencies, which is referred to as a frequency decorrelation. Attempts to detect the frequency decorrelation through a power spectra analysis of the multi-frequency Planck data have only yielded upper limits (Planck Collaboration Int. L 2017; Planck Collaboration XI 2020). This work improves on previous studies by using a new upgraded version of Planck maps. The polarization maps at frequencies 100–353 GHz used in this paper are taken from the SRo112.0 version of the Planck data processing (Delouis et al. 2019), which correct the data from systematics that plague the PR3 release. This work also improves the sensitivity to frequency decorrelation by extending the analysis from the high Galactic latitude sky regions best suited for CMB observations to brighter regions near the Galactic plane.

The paper is organized as follows. Section 2 presents the Planck data we use in this work. Section 3 introduces the reference models we used for data simulations and analysis. In Sect. 4, we determine the dust mean SED for polarization from 100 to 353 GHz. We quantify spatial variations of the dust polarization SED in Sect. 5 and the contribution of polarization angles in Sect. 6. The frequency dependence is analyzed in Sect. 7. Our results are summarized in Sect. 8.

## 2. Planck polarization data and masks

The Planck satellite observes the sky in total intensity (also referred as temperature) in the range of frequency of the electromagnetic spectrum from 30 to 857 GHz, and in polarization from 30 to 353 GHz. Data were obtained from two instruments on board the satellite: the Low Frequency Instrument (LFI, Mennella et al. 2011) and the High Frequency Instrument (HFI, Planck HFI Core Team 2011), with the Planck HFI measuring the linear polarization at 100, 143, 217, and 353 GHz (Rosset et al. 2010).

The Planck scanning strategy sampled almost all the sky pixels every six months, with alternating scan directions in successive six-month periods. The Planck mission includes five surveys, each covering a large fraction of the sky (hereafter $f_{sky}$). Maps are produced for the full-mission data set together with the survey, year, and half-mission maps, as reported in Planck Collaboration VIII (2016).

### 2.1. SRo112 maps

In this work, we use sky maps and end-to-end data simulations produced by the SRo112 software.<fnref n="1" /> This latter has been developed to improve subtraction of systematic effects. The dominant systematic effect for the polarized signal at 353 GHz in the PR3 maps is related to the poor measurement of the time transfer function of the detectors, while at lower frequencies, it is dominated by the non-linearity of the analog-to-digital converters. Both systematics have been greatly improved in a consistent way with all other known effects for the SRo112 data-set (Delouis et al. 2019).

Hereafter, we deem $Q_P(v)$ and $U_P(v)$ as the Planck polarization maps at the frequency, $v$, including the CMB, dust, and synchrotron as well as noise and systematics. We used the Planck HFI maps at 100, 143, 217, and 353 GHz. We limit our data analysis to polarization at low multipoles ($4 < l < l_{max} = 32$) and work with HEALPix pixelization (Górski et al. 2005) at $N_{side} = 32$ (i.e., map pixel size of 1.8°). In order to obtain these maps, we follow Planck Collaboration Int. XLVI (2016) and Planck Collaboration V (2020), degrading the full-resolution maps first to $N_{side} = 1024$, to ease the computation, and alongside $N_{side} = 32$, we apply the following cosine filter in harmonic space (Benabed et al. 2009):

$$f(l) = \begin{cases} 1, & l \le N_{side}; \\ \frac{1}{2}\left(1 + \sin\left(\frac{\pi}{2} \frac{l}{N_{side}}\right)\right), & N_{side} < l < 3N_{side}; \\ 0, & l \ge 3N_{side}. \end{cases} \quad (1)$$

All the maps used in this study and presented in the paper have been degraded to $N_{side} = 32$ by following this approach.

### 2.2. Subtraction of synchrotron emission

The 100 and 143 GHz maps include non-negligible synchrotron emission that we subtract to focus our data analysis on the dust emission. Hence, we define the following maps:

$$\begin{aligned} Q'_{P}(v) &= Q_{P}(v) - Q_{s}(v), \\ U'_{P}(v) &= U_{P}(v) - U_{s}(v), \end{aligned} \quad (2)$$

where $[Q_s, U_s](v)$ are estimates of the synchrotron Stokes parameters. We note that the $Q'_P(v)$ and $U'_P(v)$ maps include the CMB. Hereafter, the prime superscript is used to indicate Stokes maps where the synchrotron emission is subtracted or absent.

We used synchrotron template maps at $v_s = 30$ GHz, as obtained by the Commander component separation (Planck Collaboration X 2016) applied to PR3 Planck maps. To extrapolate from 30 to 100 and 143 GHz, we used a single spectral index, $\beta_s$, uniform over the sky, to extrapolate synchrotron polarization from 30 to $v = [100, 143]$ GHz.

$$[Q, U]_s(v) = [Q, U]_s(v_s) \cdot C_v^{CC_s} \cdot C_v^{UC_{RJ}} \cdot \left(\frac{v}{30 \text{ GHz}}\right)^{\beta_s}, \quad (3)$$

<footnote n="1">¹ The SRo112 maps are available here and the simulations here.</footnote>
<table><caption>**Table 1.** Unit conversions $C_v^{UC_{RJ},UC_K}$ and color corrections, $C_v^{CC_{d,s}}$, taken from Planck Collaboration IX (2014).</caption><thead><tr><th>ν [GHz]</th><th>100</th><th>143</th><th>217</th><th>353</th></tr></thead><tbody><tr><th>C<sub>v</sub><sup>CC<sub>d</sub></sup></th><td>1.09</td><td>1.02</td><td>1.12</td><td>1.11</td></tr><tr><th>C<sub>v</sub><sup>UC<sub>K</sub></sup></th><td>244.1</td><td>371.7</td><td>483.7</td><td>287.4</td></tr><tr><th>C<sub>v</sub><sup>UC<sub>RJ</sub></sup></th><td>1.26</td><td>1.69</td><td></td><td></td></tr><tr><th>C<sub>v</sub><sup>CC<sub>s</sub></sup></th><td>0.9797</td><td>0.9504</td><td></td><td></td></tr><tr><th>ρ<sub>v</sub></th><td>1.005</td><td>0.98</td><td>1.015</td><td>1</td></tr></tbody></table>

**Notes.** The factor $\rho(\nu)$ is the correction to the polarization efficiencies (Planck Collaboration III 2020; Planck Collaboration XI 2020).

where the $[Q, U]_s(\nu)$ maps are in $K_{CMB}$ and $[Q, U]_s(\nu_s)$ in $K_{RJ}$; $C_v^{UC_{RJ}}$ is the conversion factor from $K_{RJ}$ to $K_{CMB}$, and $C_v^{CC_s}$ is the color correction, at the frequency, $\nu$. The conversion factor values at 100 and 143 GHz are listed in Table 1.

To determine $\beta_s$, we use spectral indices derived by Martire et al. (2022), from a detailed analysis combining Planck and WMAP (Bennett et al. 2013) data at 30 and 23 GHz, respectively.

Table 4 in Martire et al. (2022) lists the spectral indices derived from EE and BB power spectra. We combine the two pairs of EE and BB spectral indices, for the two largest sky areas with $f_{sky} = 94$ and 70%, to compute a mean value, weighted by inverse squared uncertainties, $\beta_s = -3.19 \pm 0.07$.

This mean value is a reasonable approximation considering recent studies that obtained a synchrotron spectral index in polarization $\beta_s$ of: i) $-3.22 \pm 0.08$ over the southern sky (Krachmalnicoff et al. 2018); ii) $-3.17 \pm 0.06$ in the North Polar Spur (Svalheim et al. 2023); and iii) $-3.25 \pm 0.06$ within the BICEP2/Keck survey footprint (Weiland et al. 2022). The study of de la Hoz (2022) hints at evidence of spatial variations of $\beta_s$, which would need to be accounted for in future developments.

## 2.3. Sky masks and power spectra

We computed the power spectra for three sky areas presented in Fig. 1 with sky fractions, $f_{sky}$, of 80%, 90%, and 97%. We use larger values for $f_{sky}$ (than usually applied to study Galactic foregrounds) in order to increase the S/N. To mask areas of bright dust emission, we used the dust optical depth map estimated at 353 GHz by fitting a modified blackbody (MBB) spectral model to the GNILC (generalized needlet internal linear combination (GNILC, Planck Collaboration Int. XLVIII 2016) dust maps at 353, 545, 857, and the IRAS 3000 GHz map (Fixsen et al. 1999). We smoothed the map with a 5° beam before reducing the HEALPix resolution to $N_{side} = 32$. Next, we sorted the pixels by increasing amplitude to define the appropriate masks to obtain $f_{sky} = 80\%, 90\%,$ and 97%. Finally, we again smoothed the mask maps to a 5° beam resolution to apodize them and avoid edge effects on the Galactic cut contours.

The power spectra were computed by using the PolSpice estimator, which corrects for multipole-to-multipole coupling and for the mixing of the E- and B- modes due to the sky masking (Chon et al. 2004). We systematically computed the cross-power spectra to ensure we have no bias from any data noise. To compute power spectra at one given frequency, we used the so-called “half-mission” maps (hereafter “HM”). We note that most of the residual instrumental systematic effects evolve with time and are decorrelated between the two half-mission data sets.

Throughout the paper, we use $\mathcal{D}_\ell \equiv \ell(\ell+1)C_\ell/2\pi$, where $C_\ell$ is the original angular power spectrum.

# 3. Reference models and data simulations

In our analysis of the Planck data, we make use of two reference models of the dust emission and of data simulations, introduced in this section.

## 3.1. Reference models

The Planck data analysis has shown that the MBB emission law fits well the SED of the dust emission for total intensity (Planck Collaboration XI 2014; Planck Collaboration X 2016) and polarization (Planck Collaboration XI 2020). This provides a convenient and commonly used parametrization of the dust SED:

$$I_d(\nu) = C_v^{UC_K} \cdot C_v^{CC_d}(\beta_d, T_d) \cdot \rho_v \cdot \tau_{v0} \left(\frac{\nu}{\nu_0}\right)^{\beta_d} \cdot B_v(T_d) \quad (4)$$

where $B_v(T_d)$ is the Planck function; $T_d$, $\beta_d$, and $\tau_{v0}$ are the dust temperature, spectral index, and optical depth maps at the reference frequency, $\nu_0$, respectively. The MBB emission is expressed in MJy sr$^{-1}$, whereas the data are given in thermodynamic units. The conversion between the two is accomplished by two factors. The first, $C_v^{UC_K}$, is a unit conversion from MJy sr$^{-1}$ to $K_{CMB}$ for the reference spectral dependence: constant product $\nu I_v$ over the bandpass. The second, $C_v^{CC_d}(\beta_d, T_d)$, is the color correction accounting for the difference between the reference spectral dependence and the MBB spectrum. This correction depends on the MBB parameters: $\beta_d$ and $T_d$, which are taken to be equal to 1.53 and 19.6 K, respectively (Planck Collaboration XI 2020). The estimated values are presented in Table 1.

In Planck Collaboration III (2020), the polarization efficiencies have been adjusted, within the uncertainties of the ground calibration, in order to match the cosmological parameters derived from the CMB polarization with the ones obtained from CMB temperature (Planck Collaboration III 2020). The factor $\rho_v$ in Eq. (4) represents this correction to the polarization efficiencies. We use the same values as Planck Collaboration XI (2020), listed in Table 1. The uncertainty on $\rho_v$ is estimated to be 0.5% at $\nu = [100, 143, 217]$ GHz. The correction could not be estimated with the required accuracy at 353 GHz.

We applied Eq. (4) to two sets of MBB parameters derived from Planck component separations methods: i) the GNILC method using maps corrected for anisotropies of the cosmic infrared background and ii) the standard Bayesian analysis framework, implemented in the Commander code (Planck Collaboration X 2016). For GNILC, $T_d$, $\beta_d$, and $\tau_{v0}$ have been obtained by fitting a MBB model on the dust total intensity at the Planck frequencies 353, 545, and 857 GHz and IRAS 3000 GHz, while the Commander fit includes all Planck frequencies, together with the nine-year WMAP observations between 23 and 94 GHz (Bennett et al. 2013) and a 408 MHz survey map (Haslam et al. 1982).

To build our reference models, we assume that the dust SED is the same for the total intensity and polarization, an hypothesis supported by the close match between the total intensity and polarization SEDs (Planck Collaboration Int. XXII 2015; Planck Collaboration XI 2020), which we aim to further test in this paper. For GNILC, we also assume that the MBB fitted over the far-IR may be extrapolated to microwave frequencies (Planck Collaboration Int. XVII 2014). Within this framework, the Stokes parameters $Q_d(\nu)$ and $U_d(\nu)$ at $\nu =$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Apodized masks for $f_{sky} = 80, 90$, and $97\%$ (from left to right).</figcaption></figure>

[100, 143, 217] GHz may be computed from the Planck maps at $v_0 = 353$ GHz:

$$ 
\begin{aligned}
Q_d(v) &= \frac{I_d(v)}{I_d(v_0)} \cdot (Q_P(v_0) - Q_{P,CMB}), \\
U_d(v) &= \frac{I_d(v)}{I_d(v_0)} \cdot (U_P(v_0) - U_{P,CMB}),
\end{aligned} 
\qquad (5) $$

where $Q_{P,CMB}$, and $U_{P,CMB}$ are CMB Stokes maps from Planck. We use the Planck SMICA CMB maps (Planck Collaboration IV 2020). Our Commander and GNILC models are low-resolution versions of dust models in PySM (Thorne et al. 2017; Zonca et al. 2021), which are commonly used by the CMB community.

## 3.2. Data simulations

In order to account for the Planck satellite noise, including instrumental systematics and the CMB signal, we computed a set of simulated maps as:

$$ 
\begin{aligned}
Q'_{\text{sim}}(v) &= Q_d(v) + Q_{\text{noise+syst}}(v) + Q_{\text{CMB}}, \\
U'_{\text{sim}}(v) &= U_d(v) + U_{\text{noise+syst}}(v) + U_{\text{CMB}},
\end{aligned} 
\qquad (6) $$

where $[Q_{\text{noise+syst}}, U_{\text{noise+syst}}]$ are the 200 SR0112 simulations to which we subtract the Planck sky model that contain synchrotron, dust, and CMB in order to isolate Planck noise and systematics, (see Delouis et al. 2019, for more details). In addition, $[Q_{\text{CMB}}, U_{\text{CMB}}]$ are 200 independent realizations of the CMB computed from the theoretical power spectra of the best-fit ΛCDM model to the Planck data<fnref n="2" /> (Planck Collaboration VI 2020), with the r parameter equal to 0. Thus, Eq. (6) yields 200 realizations of the reference model maps.

In Fig. 2, we present one realization of the Commander $Q'_{\text{sim}}$ and $U'_{\text{sim}}$ maps at 100, 143, and 217 GHz. In Fig. 3, the power spectra of the dust model Commander, the synchrotron template at the two lowest frequencies, the CMB, and a simulation of the SR0112 noise plus systematics are compared. The mean values $\langle D_t \rangle_{t=[4,32]}$ for EE (diamonds) and BB (squares) spectra are plotted versus frequency for the three values of $f_{sky}$.

## 4. Dust mean SED in polarization

In this section, we derive the mean SED of dust polarization. We use Planck polarization maps at 100, 143, 217, and 353 GHz. The SED values are normalized to the reference frequency, $v_0 = 353$ GHz.

<figure><img src="image_3.png" /><figcaption>Fig. 2. $Q'_{\text{sim}}$ and $U'_{\text{sim}}$ reference maps for 100 GHz (top), 143 GHz (middle), and 217 GHz (bottom) as obtained for the Commander reference model.</figcaption></figure>

We follow earlier studies (Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2020) using cross power spectra to determine the dust mean polarized SED for Planck data $\gamma_p(v)$, the reference models $\gamma_d(v)$, and the simulations $\gamma_{\text{sim}}(v)$ normalized to $v_0$ as:

$$ 
\begin{aligned}
\gamma_{P}^{\text{XX}}(v) &= \rho_{v} \left\langle \frac{\mathcal{D}_{\ell}^{\text{XX}}(v \times v_0) - \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}}}{\mathcal{D}_{\ell}^{\text{XX}}(v_0 \times v_0) - \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}}} \right\rangle_{\ell_{\text{min}}, \ell_{\text{max}}}, \\
\gamma_{d}^{\text{XX}}(v) &= \left\langle \frac{\mathcal{D}_{\ell, \text{d}}^{\text{XX}}(v \times v_0)}{\mathcal{D}_{\ell, \text{d}}^{\text{XX}}(v_0 \times v_0)} \right\rangle_{\ell_{\text{min}}, \ell_{\text{max}}}, \\
\gamma_{\text{sim}}^{\text{XX}}(v) &= \rho_{v} \left\langle \frac{\mathcal{D}_{\ell, \text{sim}}^{\text{XX}}(v \times v_0) - \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}}}{\mathcal{D}_{\ell, \text{sim}}^{\text{XX}}(v_0 \times v_0) - \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}}} \right\rangle_{\ell_{\text{min}}, \ell_{\text{max}}},
\end{aligned} 
\qquad (7) $$

where

$$ 
\begin{aligned}
\mathcal{D}_{\ell}^{\text{XX}}(v \times v_0) &= [Q'_{P}^{\text{HM1}}(v), U'_{P}^{\text{HM1}}(v)] \times [Q'_{P}^{\text{HM2}}(v_0), U'_{P}^{\text{HM2}}(v_0)], \\
\mathcal{D}_{\ell, \text{d}}^{\text{XX}}(v \times v_0) &= [Q_{\text{d}}^{\text{HM1}}(v), U_{\text{d}}^{\text{HM1}}(v)] \times [Q_{\text{d}}^{\text{HM2}}(v_0), U_{\text{d}}^{\text{HM2}}(v_0)], \\
\mathcal{D}_{\ell, \text{sim}}^{\text{XX}}(v \times v_0) &= [Q'_{\text{sim}}^{\text{HM1}}(v), U'_{\text{sim}}^{\text{HM1}}(v)] \times [Q'_{\text{sim}}^{\text{HM2}}(v_0), U'_{\text{sim}}^{\text{HM2}}(v_0)],
\end{aligned} 
\qquad (8) $$

<footnote n="2">² We used the data file COM_PowerSpect_CMB-base-plikHM-TTTEE-Lowl-lowE-lensing-minimum-theory_R3.01.txt available on the Planck Legacy archive.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Amplitudes, averaged for $4 < l < 32$, of the $D_l^{EE}$ (diamond) and $BB$ (square) power spectra versus frequency. Each plot presents the spectra of the Commander dust model (gray), our synchrotron estimate at 100 and 143 GHz (black), the CMB (orange), and the SRo112 noise plus systematics (light brown). The three plots from left to right correspond to $f_{sky} = 80, 90$, and $97\%$.</figcaption></figure>

<table><caption>**Table 2.** Dust mean SED values for the Planck polarization data $\gamma_P(\nu)$ and the reference models $\gamma_d(\nu)$.</caption><thead><tr><th>f<sub>sky</sub></th><th colspan="3">80%</th><th colspan="3">90%</th><th colspan="3">97%</th></tr><tr><th>ν[GHz]</th><th>100</th><th>143</th><th>217</th><th>100</th><th>143</th><th>217</th><th>100</th><th>143</th><th>217</th></tr></thead><tbody><tr><th>γ<sup>EE</sup><sub>P</sub></th><td>0.0181</td><td>0.0400</td><td>0.1276</td><td>0.0178</td><td>0.0392</td><td>0.1274</td><td>0.0188</td><td>0.0411</td><td>0.1314</td></tr><tr><th>γ<sup>EE</sup><sub>d</sub>Commander</th><td>0.0187</td><td>0.0388</td><td>0.1291</td><td>0.0185</td><td>0.0382</td><td>0.1286</td><td>0.0194</td><td>0.0395</td><td>0.1315</td></tr><tr><th>γ<sup>EE</sup><sub>d</sub>GNILC</th><td>0.0177</td><td>0.0373</td><td>0.1265</td><td>0.0182</td><td>0.0379</td><td>0.1279</td><td>0.0209</td><td>0.0416</td><td>0.1347</td></tr><tr><th>error</th><td>0.0005</td><td>0.0003</td><td>0.0005</td><td>0.0002</td><td>0.0002</td><td>0.0003</td><td>0.0002</td><td>0.0001</td><td>0.0002</td></tr><tr><th>γ<sup>BB</sup><sub>P</sub></th><td>0.0180</td><td>0.0394</td><td>0.1269</td><td>0.0173</td><td>0.0391</td><td>0.1269</td><td>0.0180</td><td>0.0393</td><td>0.1276</td></tr><tr><th>γ<sup>BB</sup><sub>d</sub>Commander</th><td>0.0186</td><td>0.0387</td><td>0.1291</td><td>0.0182</td><td>0.0382</td><td>0.1285</td><td>0.0183</td><td>0.0384</td><td>0.1292</td></tr><tr><th>γ<sup>BB</sup><sub>d</sub>GNILC</th><td>0.0174</td><td>0.0371</td><td>0.1262</td><td>0.0181</td><td>0.0379</td><td>0.1279</td><td>0.0187</td><td>0.0388</td><td>0.1297</td></tr><tr><th>error</th><td>0.0005</td><td>0.0005</td><td>0.0008</td><td>0.0003</td><td>0.0003</td><td>0.0004</td><td>0.0002</td><td>0.0002</td><td>0.0003</td></tr></tbody></table>

**Notes.** The uncertainty is estimated as standard deviation over 200 simulations. For 100 and 143 GHz the uncertainty accounts for the error associated with the synchrotron template subtraction.

and $D_{l,CMB}^{XX}$ is the CMB power spectrum for the $\Lambda$CDM fiducial Planck model (Planck Collaboration VI 2020) with $XX \in [EE, BB]$. $D_{l}^{XX}(v_0 \times v_0)$ is computed with $v = v_0$ in Eq. (8). The symbol $\times$ indicates the cross-power spectrum operator and $\langle\rangle$ the arithmetic mean over the $l$-range from $l_{min} = 4$ to $l_{max} = 32$. We note that in the denominator of the first and third equations of Eq. (7), we omit $\rho_{v_0}$ because it is equal to 1. Table 2 lists the SED values $\gamma_P(\nu)$ for the Planck data and $\gamma_d(\nu)$ for the reference models, for both EE and BB and for the three sky areas. The uncertainties $\sigma_{\gamma_P}(\nu)$ are derived from the standard deviation of the 200 Planck simulations. The uncertainty on the polarization efficiencies, $\rho_v$, and the synchrotron spectral index, $\beta_s$, are propagated through our analysis and their contribution to the total error-bar are added to $\sigma_{\gamma_P}$.

The dust SEDs Planck $\gamma_P$ are presented in Fig. 4. The values are normalized to a MBB SED with $\beta_d = 1.53$ and $T_d = 19.6$ K (Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XXII 2015). We note that the value at 353 GHz is always equal to unity because it is the refer-

ence frequency considered for our analysis, see Eq. (7). The color corrections $C_v^{CC_d}(\beta_d = 1.53$ and $T_d = 19.6$ K) are those listed in Table 1. The figure shows that the $\gamma_P(\nu)$ values are consistent with the MBB used for normalization within 5%, in agreement with previous results (Planck Collaboration Int. XXII 2015; Planck Collaboration XI 2020). Some of the differences could be due to systematics that plague PR3 data. We note that subtracting the synchrotron emission primarily affects the value of the EE signal to 100 GHz. Although its level is low, around 2–3%. For the mask at 97%, we observe a significant difference between $\gamma_P^{EE}(\nu)$ and $\gamma_P^{BB}(\nu)$ that we interpret as variations of the SED within the beam, which do not average in the same way for EE and BB. These variations depend also on the mask and could prove more important when the analysis includes a significant part of the Galactic plane.

Figure 5 shows the ratio between $\gamma_d$ for the Commander (top) and GNILC (bottom) models and $\gamma_P$. For Commander, we observe a very close match between $\gamma_P(\nu)$ and $\gamma_d(\nu)$ at 143 and 217 GHz for $f_{sky} = 80$ and 90% and both EE and BB signals.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Dust-mean SED $\gamma_P(\nu)$ normalized by a modified black body function $I_d(\nu)$ with fixed $\beta_d = 1.53$ and $T_d = 19.6$ K as given by Planck Collaboration XI (2020) and accounting for color corrections. The EE and BB values are shown in red and blue, respectively. From left to right: Results for $f_{sky} = [80, 90, 97]\%$.</figcaption></figure>

For GNILC, the match is not as good – we observe in particular a significant difference for $f_{sky} = 80\%$ at 100 and 143 GHz. For $f_{sky} = 97\%$, we observe a difference between the EE and BB values at 100 GHz. Apart for this latter difference, we observe that $\gamma_d(\nu)$ values are consistent within 5% with the $\gamma_P(\nu)$ values. All the differences between $\gamma_P$ and $\gamma_d$ may be due to averaging effects along the line of sight.

This analysis shows that assuming a reference model based on total intensity data could not completely reproduce the SED observed in polarization data, thus biasing any CMB polarization E-modes and B-modes signals if used as reference for component separation methods. In order to ensure an unbiased detection of the CMB polarization at the precision required from future CMB experiments, we need to address these very small spatial SED variations in the dust polarization emission.

# 5. Spatial variations in the polarization SED

In this section, we characterize the spatial variations in the dust SED in polarization and quantify the degree of correlation with variations in the dust SED in total intensity.

## 5.1. Residual maps

To quantify the spatial variations in the dust SED, we compute the differences between the Planck frequency maps at 100, 143, and 217 GHz and the 353 GHz scaled by the mean dust SED $\gamma_P(\nu)$ from Sect. 4, which we approximate as the value of $\gamma_P^{EE}(\nu)$ because the dust polarization is dominated by E-modes. Figure 4 shows that $\gamma_P^{EE}(\nu) \approx \gamma_P^{BB}(\nu)$ for $f_{sky} = 80\%$, 90%. A difference between the two coefficients $\gamma_P(\nu)$ is instead detected for $f_{sky} = 97\%$. This difference indicates a correlation between dust emission properties and the structure of the magnetized interstellar medium. However, we have checked that our approximation also does not significantly impact the results for this $f_{sky}$. The three maps, $R_Q(\nu)$ and $R_U(\nu)$, are then computed as:

$$R_Q(\nu) = Q'_P(\nu) - \gamma_P(\nu) \cdot Q_P(\nu_0), \qquad (9)$$

$$R_U(\nu) = U'_P(\nu) - \gamma_P(\nu) \cdot U_P(\nu_0),$$

with $\nu = 100, 143$ and $217$ GHz, and $\nu_0 = 353$ GHz, which we hereafter refer to as the residual maps and shown in Fig. 6. For comparison, we also compute difference maps for the reference models (see Eq. (10)) and simulations maps replacing $\gamma_P(\nu)$ by $\gamma_d(\nu)$ in Eq. (9):

$$R_{Q_d}(\nu) = Q_d(\nu) - \gamma_d(\nu) \cdot Q_d(\nu_0), \\ R_{U_d}(\nu) = U_d(\nu) - \gamma_d(\nu) \cdot U_d(\nu_0), \qquad (10)$$

$$R_{Q'_{sim}}(\nu) = Q'_{sim}(\nu) - \gamma_d(\nu) \cdot Q'_{sim}(\nu_0),$$

$$R_{U'_{sim}}(\nu) = U'_{sim}(\nu) - \gamma_d(\nu) \cdot U'_{sim}(\nu_0).$$

We note that for the purposes of illustration, all the residual maps shown in Figs. 6–8 consider a $\gamma_P(\nu)$ and $\gamma_d(\nu)$ in Eq. (9) computed for $f_{sky} = 90\%$.

One set of simulation maps [$R_{Q'_{sim}}(\nu)$, $R_{U'_{sim}}(\nu)$] for the Commander model is displayed in Fig. 7. The comparison between Figs. 6 and 7 is hampered by Planck data noise but we may notice some common features and some differences. In Fig. 6, we also note differences between frequencies among the $R_Q$ and $R_U$ residual maps, in contrast to what is observed for $R_{Q'_{sim}}$ and $R_{U'_{sim}}$ in Fig. 7. To illustrate that some of the structures are above data noise and uncorrected systematics, in Fig. 8 we present the independent $R_Q$ and $R_U$ residual half-mission HM1 and HM2 maps estimated at 217 GHz. This qualitative examination of the maps suggests that we observe SED variations in the residual $R_Q(\nu)$, $R_U(\nu)$ maps, which cannot be simply attributed to SED variations in the total intensity, as assumed in the reference models.

## 5.2. Power-spectra analysis

To characterize the SED variations, we computed power spectra of residual maps and cross-power spectra with the reference models. We computed two cross-power spectra between (i) independent residual maps [$R_Q^{HM1, HM2}(\nu)$, $R_U^{HM1, HM2}(\nu)$] and (ii) between [$R_Q^{HM1}(\nu)$, $R_U^{HM1}(\nu)$] and the reference model residuals [$R_{Qd}^{HM2}(\nu)$, $R_{Ud}^{HM2}(\nu)$]. The 200 simulations are used to assess the impact of data noise plus uncorrected systematics and of the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Dust mean SED $\gamma_d(\nu)$ computed for the reference models: Commander (top row) and GNILC (bottom row). EE and BB values are shown in red and blue, respectively. From left to right: Results for $f_{sky} = [80, 90, 97]\%$. All values are normalized by the $\gamma_p(\nu)$ values obtained from the data and shown in Fig. 4.</figcaption></figure>

CMB on our analysis.

$$ 
\begin{aligned}
\mathcal{D}_{\ell, \text{res}}^{\text{XX}}(\nu) &= [R_Q^{\text{HM1}}(\nu), R_U^{\text{HM1}}(\nu)] \times [R_Q^{\text{HM2}}(\nu), R_U^{\text{HM2}}(\nu)], \\
\mathcal{D}_{\ell, \text{d, res}}^{\text{XX}}(\nu) &= [R_Q^{\text{HM1}}(\nu), R_U^{\text{HM1}}(\nu)] \times [R_Q^{\text{HM2}}(\nu), R_U^{\text{HM2}}(\nu)], \\
\mathcal{D}_{\ell, \text{sim, res}}^{\text{XX}}(\nu) &= [R_{Q'_{\text{sim}}}^{\text{HM1}}(\nu), R_{U'_{\text{sim}}}^{\text{HM1}}(\nu)] \times [R_{Q'_{\text{sim}}}^{\text{HM2}}(\nu), R_{U'_{\text{sim}}}^{\text{HM2}}(\nu)],
\end{aligned} 
\quad (11) $$

where $XX \in [EE, BB]$. Figure 9 shows both sets of cross power spectra for EE (red) and BB (blue) at 100, 143, and 217 GHz for $f_{sky} = 90\%$. The data points over the $\ell$-range 4–32 are compared with the $1\sigma$ and $2\sigma$ dispersion (dark and light gray shades) obtained with the 200 simulations, including noise and CMB anisotropies (see Sect. 3.2). In all plots, individual data points have a low S/N and it is necessary to average them to quantify mean amplitudes and their frequency dependence.

The averaged cross power spectra values are computed as:

$$ 
\begin{aligned}
A_{\text{res}}(\nu) &= \langle \mathcal{D}_{\ell, \text{res}}^{\text{XX}}(\nu) - (1 - \gamma_{\text{P}}(\nu))^2 \cdot \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}} \rangle_{\ell_{\text{min}}, \ell_{\text{max}}}, \\
A_{\text{res,d}}(\nu) &= \langle \mathcal{D}_{\ell, \text{d, res}}^{\text{XX}}(\nu) \rangle_{\ell_{\text{min}}, \ell_{\text{max}}}, \\
A_{\text{res,sim}}(\nu) &= \langle \mathcal{D}_{\ell, \text{res,sim}}^{\text{XX}}(\nu) - (1 - \gamma_{\text{P}}(\nu))^2 \cdot \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}} \rangle_{\ell_{\text{min}}, \ell_{\text{max}}},
\end{aligned} 
\quad (12) $$

with $XX \in [EE, BB]$. Notice that the CMB has been partially subtracted in Eqs. (9) and (10). The remaining contribution is accounted for by the term $(1 - \gamma_{\text{P}}(\nu))^2 \cdot \mathcal{D}_{\ell, \text{CMB}}^{\text{XX}}$. In Fig. 10, we plot the amplitudes $A_{\text{res}}(\nu)$ normalized by the mean polarized intensity at reference frequency over each sky area defined as:

$$ P_{v_0} = \sqrt{|Q_{HM1} \cdot Q_{HM2} + U_{HM1} \cdot U_{HM2}|}, \quad (13) $$

where HM1 and HM2 refer to the two half-mission maps. We checked that the bias on this estimator of the polarized intensity is negligible at 353 GHz for $N_{\text{side}} = 32$. This normalization allows us to compare the results for the three different masks.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Residual maps $R_Q$ (left) and $R_U$ (right) at 100, 143, and 217 GHz (shown from top to bottom).</figcaption></figure>

Absolute values of the residuals may be obtained by scaling the figure data points by $\langle P_{v_0}^2 \rangle = [2453.6, 5950.5, 13814.3]$ K$_{\text{CMB}}$ for $f_{\text{sky}} = 80\%, 90\%$, and $97\%$ respectively. We discuss the frequency dependence later, here we focus on the comparison with the models. Figure 11 shows the dust models residuals $A_{\text{res,d}}(\nu)$ normalized by $A_{\text{res}}(\nu)$. The mean amplitudes $A_{\text{res,sim}}(\nu)$ are also computed for each of the 200 simulations and their dispersion provides the error bars. For the 80 and 90% masks, the $A_{\text{res}}(\nu)$
<figure><img src="image_8.png" /><figcaption>Fig. 7. Residual maps for one realization of $R_{Q'_{sim}}$ (left) and $R_{U'_{sim}}$ (right) obtained for the Commander reference model, at 100, 143, 217 GHz respectively (shown from top to bottom).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Residual maps $R_Q$ (left) and $R_U$ (right) at 217 GHz for two data sets HM1 (top), HM2 (bottom), respectively (shown from top to bottom).</figcaption></figure>

are larger by a factor of about ten larger than $A_{res,d}(v)$. For these, uncertainties associated with the chance correlation between the CMB and the dust residuals may thus be somewhat underestimated. For the 97% mask, this difference is smaller. At 143 and 217 GHz, the amplitudes even coincide for the GNILC model. These results show that the reference model only accounts for a minor fraction of the total polarization SED variations in the 80 and 90% masks but a much more significant one for the brighter dust emission in the Galactic plane. Furthermore, Fig. 10 highlights the significant variation of BB residuals at 100 GHz toward higher galactic latitudes. Although both the EE and BB variations are very small, in the case of $f_{sky} = 80\%$, we have $A_{res}^{BB}$ representing 0.4% of the total amplitude with regard to the 0.1% detected for $A_{res}^{EE}$. This feature does not match the residuals extrapolated from the models considered (see Fig. 11 for comparison).

To give a more quantitative estimate, we translate the amplitude, $A_{res}(v)$, in terms of an effective dispersion of the dust spectral index, $\sigma_{\beta}$, using the following formula based on a first-order expansion of the MBB emission law in Mangilli et al. (2021):

$$ \sigma_{\beta}(v) = (A_{\text{res}}^{\text{EE}}(v) + A_{\text{res}}^{\text{BB}}(v))^{0.5} \cdot \left(\gamma_{\text{P}}(v) \cdot P_{v0} \cdot \left|\ln\left(\frac{v}{v_0}\right)\right|\right)^{-1} . \quad (14) $$

The values of $\sigma_{\beta}$ are listed in Table 3 for the three sky areas and frequencies. These numbers (of about 0.1) must be treated as an estimate of the variance of residuals in terms of pure variation of $\beta$. These results are consistent with previous studies at low latitude (Planck Collaboration Int. XLVIII 2016).

## 6. Residuals from polarization angles

In this section, we quantify the variations in the dust polarization angles as a function of frequency. To do this we introduce the

following maps:

$$ \tilde{Q}(v) = \gamma_{\text{P}}(v) \cdot P(v_0) \times \cos(2\psi(v)), \quad (15) $$

$$ \tilde{U}(v) = \gamma_{\text{P}}(v) \cdot P(v_0) \times \sin(2\psi(v)), $$

where $P(v_0)$ is the polarized intensity as estimated in Eq. (13), and $\psi(v) = \frac{1}{2} \arctan(\frac{U}{Q})$ is the polarization angle. To assess the contribution of variations in terms of $\psi$ to the total $A_{res}(v)$ values, we also computed residual maps for $\tilde{Q}(v)$ and $\tilde{U}(v)$ as:

$$ \tilde{R}_{Q}(v) = \tilde{Q}(v) - \gamma_{\text{P}}(v) \cdot \tilde{Q}(v_0), \quad (16) $$

$$ \tilde{R}_{U}(v) = \tilde{U}(v) - \gamma_{\text{P}}(v) \cdot \tilde{U}(v_0). $$

The power spectra analysis is performed with these residual maps as described in the previous section:

$$ \tilde{D}_{\ell, \text{res}}^{\text{XX}}(v) = [\tilde{R}_{Q}^{\text{HM1}}(v), \tilde{R}_{U}^{\text{HM1}}(v)] \times [\tilde{R}_{Q}^{\text{HM2}}(v), \tilde{R}_{U}^{\text{HM2}}(v)], \quad (17) $$

$$ \tilde{D}_{\ell, \text{sim, res}}^{\text{XX}}(v) = [\tilde{R}_{Q'_{\text{sim}}}^{\text{HM1}}(v), \tilde{R}_{U'_{\text{sim}}}^{\text{HM1}}(v)] \times [\tilde{R}_{Q'_{\text{sim}}}^{\text{HM2}}(v), \tilde{R}_{U'_{\text{sim}}}^{\text{HM2}}(v)]. $$

And then the averaged amplitudes are defined as:

$$ \tilde{A}_{\text{res}}(v) = \langle \tilde{D}_{\ell, \text{res}}^{\text{XX}}(v) - \langle \tilde{D}_{\ell, \text{sim, res}}^{\text{XX}}(v) \rangle_{\text{sim}} \rangle_{\ell_{\text{min}}, \ell_{\text{max}}}, \quad (18) $$

where the second term is used for the CMB and residual noise and systematics debiasing. As in Sect. 5.2 for $A_{res}$, the uncertainties on $\tilde{A}_{res}(v)$ are derived from the dispersion of the values obtained for the simulations.

Figure 12 compares the $\tilde{A}_{res}(v)$ values to $A_{res}(v)$. We find that variations of the polarization angle contribute significantly to the total polarization residuals. The differences observed with $f_{sky}$ may result from the integration along the line of sight or from the limit of our debiasing method (or both). Indeed, for decreasing Galactic latitudes, the line of sight crosses an increasing number of coherent turbulent cells and the impact of the magnetic field structure on observed polarization angles may average out.

To express the amplitude, $\tilde{A}_{res}(v)$, in terms of an effective angle variation, $\delta_{\psi}$, we assume that the residuals results from a systematic angle change over the full sky area. This calculation is solely indicative because we do not believe that this
<figure><img src="image_10.png" /><figcaption>Fig. 9. Cross-power spectra correlation $D_{\ell}^{XX}(\nu)$ and $D_{\ell,d}^{XX}(\nu)$, indicated as HM1 × HM2 and HM1 × Ref model HM2, respectively, in the labels. Results for $f_{sky} = 90\%$ and the Commander reference model are shown. First two rows show the $D_{\ell}^{EE}$ and the last two the $D_{\ell}^{BB}$ cross power spectra, respectively (from top to bottom). Results at 100, 143, and 217 GHz are shown from left to right. The dark and light gray shades represent the $1\sigma$ and $2\sigma$ standard deviation computed from the data simulations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Amplitudes $A_{res}(\nu)$ for EE (red) and BB (blue), respectively. These values are averages over the $\ell$-range 4–32 and normalized by the mean of $P_{y0}^2$. Results for $f_{sky} = [80, 90, 97]\%$ are shown (from left to right).</figcaption></figure>

is a valid assumption. Using equations detailed in Abitbol et al. (2016), we obtain:

$$ \sin \delta_{\psi}(\nu) = 0.5 \left( \tilde{A}_{\text{res}}^{\text{EE}}(\nu) + \tilde{A}_{\text{res}}^{\text{BB}}(\nu) \right)^{0.5} \gamma_{\text{P}}^{-1}(\nu) \langle (P_{y0}^2) \rangle^{-0.5}. \quad (19) $$

The values of $\delta_{\psi}$ are listed in Table 3 for three sky areas and frequencies. These values are small but still larger than the $1^\circ$ uncertainty on the ground calibration of the Planck absolute polarization angle (Rosset et al. 2010) at both 100 and 143 GHz. The analysis of the Planck data confirms this upper limit of $1^\circ$ (see Fig. 20 in Delouis et al. 2019).

## 7. Frequency dependence

In this section, we discuss the frequency dependence of the amplitudes of the residuals power spectra plotted in Fig. 10. We follow earlier studies (Chluba et al. 2017; Désert 2022) using a Taylor expansion of the MBB emission law to model the residuals. The moment expansion has been applied to dust power spectra in total intensity by Mangilli et al. (2021) and to dust B-modes power spectra considered as an intensity by Azzoni et al. (2021) and Vacher et al. (2022b). Within this framework, the power spectra of the residual maps are modeled
<figure><img src="image_12.png" /><figcaption>Fig. 11. Residual amplitudes $A_{res,d}(v)$ obtained for the reference models Commander (top) and GNILC (bottom). Red and blue colors represent the EE and BB results, respectively. These values are estimated by averaging between $l$-range 4–32 and normalizing to the amplitudes of the residuals obtained from the Planck data $A_{res}(v)$. See Eq. (12) for details. Results for $f_{sky} = [80, 90, 97]\%$ are shown from left to right. The absolute value of BB negative results found at 100 GHz for $f_{sky} = 80\%$ is shown in gray.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Amplitudes of the residuals $\tilde{A}_{res}(v)$ obtained by isolating the effect of the polarization angle variation are shown for EE in red and BB in blue (see Eq. (18)). These values are normalized to the residual amplitudes obtained from the Planck data analysis $A_{res}(v)$. Values for $f_{sky} = 80, 90$, and $97\%$ are shown from left to right.</figcaption></figure>

using the first-order expansion as:

$$ \begin{aligned} \mathcal{D}_t(v) = \gamma_P(v)^2 \cdot \left\{ \right. & \text{1st order } \beta \left\{ +\mathcal{D}_l^{\omega_l^\beta \times \omega_l^\beta} \ln\left(\frac{v}{v_0}\right)^2, \right. \\ & \text{1st order } T \left\{ +\mathcal{D}_l^{\omega_l^T \times \omega_l^T} (\Theta_v(T_d) - \Theta_{v_0}(T_d))^2, \right. \\ & \text{1st order } T \times \beta \left\{ +2\mathcal{D}_l^{\omega_l^\beta \times \omega_l^T} \ln\left(\frac{v}{v_0}\right) \cdot (\Theta_v(T_d) - \Theta_{v_0}(T_d)) \right\}. \end{aligned} \quad (20) $$

Here, $\mathcal{D}_l^{\alpha \times \beta}$ represents three moment coefficients introduced by Mangilli et al. (2021) and Vacher et al. (2022b), which are associated

with spatial variations of the dust spectral index and temperature and correlated variations of these two MBB parameters. The $\Theta_v(T_d)$ function is the derivative of the logarithm of the black-body spectrum with respect to temperature $T_d$:

$$ \Theta_v(T_d) = \frac{x}{T_d} \frac{e^x}{e^x - 1}, \quad (21) $$

where $x = \frac{h\nu}{k_B T_d}$. The temperature, $T_d$, is the one used to normalize $\gamma_P(v)$ in Fig. 4.

Ichiki et al. (2019) – and even more extensively Vacher et al. (2023) – have extended the moment expansion to polarization to model the frequency dependence of Stokes Q and U maps. In
<table><caption>**Table 3.** Effective dispersion of the dust spectral index $\sigma_{\beta}$ and effective angle variation $\delta_{\psi}$.</caption><thead><tr><th rowspan="2">f<sub>sky</sub></th><th colspan="3">80%</th><th colspan="3">90%</th><th colspan="3">97%</th></tr><tr><th>100</th><th>143</th><th>217</th><th>100</th><th>143</th><th>217</th><th>100</th><th>143</th><th>217</th></tr></thead><tbody><tr><th>ν[GHz]</th><td>100</td><td>143</td><td>217</td><td>100</td><td>143</td><td>217</td><td>100</td><td>143</td><td>217</td></tr><tr><th>σ<sub>β</sub></th><td>0.15±0.03</td><td>0.12 ± 0.01</td><td>0.09±0.02</td><td>0.12±0.02</td><td>0.09±0.01</td><td>0.08±0.01</td><td>0.10±0.01</td><td>0.07 ± 0.01</td><td>0.07 ± 0.01</td></tr><tr><th>δ<sub>ψ</sub></th><td>2.71±0.46</td><td>1.60 ± 0.25</td><td>0.81±0.19</td><td>2.43±0.29</td><td>1.13±0.15</td><td>0.66±0.11</td><td>2.33±0.20</td><td>0.83 ± 0.09</td><td>0.48 ± 0.07</td></tr></tbody></table>

**Notes.** A systematic uncertainty on the polarization angle of 1° must be considered as upper limit on the Planck polarization absolute accuracy (Rosset et al. 2010).

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Residual amplitudes $A_{\text{res}}(\nu)$, (dark markers) EE (red) and BB (blue) for $f_{\text{sky}}$ = [80, 90, 97]%, are shown from left to right. A fit to the data EE accounting for the three components of the moment expansion modeling separately, is described in Eq. (20) and is shown as a dashed line. A cross-power spectra correlation between pair of frequencies: 100 × 143, 100 × 217 and 143 × 217 is shown in smaller and lightened markers.</figcaption></figure>

the formalism introduced by Vacher et al. (2023), the moments are spin-2 objects that characterize the frequency dependence of both polarized intensity and polarization angle. This formalism has not yet been applied to EE and BB power spectra. Qualitatively, we would expect that variations of polarization angles induce an exchange of power between E- and B-modes that is not symmetric due to the E-B asymmetry of dust polarization. The frequency dependence of EE and BB power spectra of residuals maps are thus coupled – but not necessarily identical.

In Fig. 13, we present fits of the $A_{\text{res}}^{\text{EE}}(\nu)$ with the three terms in Eq. (20), which we refer to as $f_1$, $f_2$ and $f_3$ for $\beta$, $T$ and $T \times \beta$, respectively. Over the frequency range of our data analysis from 100 to 217 GHz, these three functions differ only slightly. The $f_1$ function provides the best fit, especially for $f_{\text{sky}}$ = 97%, but the fits do not allow us to disentangle the contributions from variations of $\beta_d$, $T_d$ and their correlation to the power spectra of the residual maps.

It is satisfactory to find that the moments' expansion to the first order provides a good model for the EE power-spectra of residual maps. This result suggests that the residuals follow mainly from variations of dust spectral parameters. However, we can see in the figure that the fits would not be as good for the BB amplitudes, in particular for $f_{\text{sky}}$ = 80%, where the amplitudes do not increase with increasing frequency.

Figure 13 includes amplitudes derived from the 100 × 143, 100 × 217, and 143 × 217 cross-spectra that are plotted with smaller symbols. These data points were not used in the fit because their values depend on frequency decorrelation. We do find that the residual maps at the three frequencies are not fully correlated. The 143 × 217 (and to a lesser extent the 100 × 217)

amplitudes lie under the model fit. Part of this mismatch could result from a miscalibration of the absolute polarization angle at 217 GHz. To investigate this possibility, we performed the following test. We rotated the reference frame of the $Q'_p$, and $U'_p$ maps at 217 GHz by an angle $\psi = \pm 1^\circ$, which is an upper limit based on the Planck polarization absolute angle uncertainty (Delouis et al. 2019). We repeated our analysis with these rotated Stokes maps at 217 GHz for both the Planck data and the simulations to correct for the impact of the rotation on the CMB signal. We found that the amplitudes of 100 × 217 and 143 × 217 do not change in any marked way for either EE and BB, nor for the three masks considered. We conclude from this test that the mismatch between the spectral model and the cross-spectra amplitudes cannot be explained by a miscalibration of the polarization angle in the Planck 217 GHz data.

## 8. Conclusions

Power spectra of the Planck data were used to characterize spatial variations of the polarized dust SED. We improved the sensitivity of previous studies by using the newly released SRo112 maps and extending the analysis to regions near the Galactic plane. Our analysis focuses on the lowest multipoles between $l = 4$ and 32, and three sky areas with $f_{\text{sky}}$ = 80%, 90%, and 97%. Maps of MBB parameters from the Commander and GNILC component separation methods applied to the Planck total intensity data are used as reference models. The main results of our analysis are as follows.

- We confirm earlier studies finding that the mean SED for dust polarization from 100 to 353 GHz is very close to that
for the total intensity, as well as to a MBB spectrum with a spectral index of $\beta_d = 1.53$ for an assumed dust temperature of $T_d = 19.6$ K.

- The mean SED and the 353 GHz *Q* and *U* maps are used to compute residual maps at 100, 143, and 217 GHz, which quantify spatial variations of the dust polarization SED. Residuals are detected at the three frequencies for the three sky areas. The *EE* and *BB* spectra of the residual maps are of comparable amplitude. They do not reproduce the *E* – *B* asymmetry that has been observed for the total dust power.

- The residual maps are correlated with the reference Commander and GNILC models, but this correlation accounts for only a fraction of the residuals amplitude. Furthermore, we find that this fraction decreases toward high Galactic latitudes (i.e., for a decreasing $f_{sky}$). This result shows that models based on total intensity data are underestimating the complexity of dust polarized CMB foreground. Possibly, future developments of Commander and GNILC models, also accounting for the latest releases of the *Planck* data, could improve the comparison with polarized data.

- To gain insights into the origin of SED variations, we quantify the variations in the polarization angle. For $f_{sky}$ = 80% and 90%, we find that the contribution of polarization angles to the residuals is dominant, in particular for the *BB* signal. These results emphasize the importance of considering the geometrical properties of Galactic polarization in component separation.

- The frequency dependence of the *EE* and *BB* residual amplitudes yields further insights. We find that the moments' expansion to the first order of the MBB spectrum provides a good fit to the *EE* amplitudes. This result suggests that the residuals follow mainly from variations of dust spectral parameters (temperature and spectral index). However, this conclusion is challenged by the *BB* results, particularly for $f_{sky}$ = 80%, and by cross-spectra that show that the residuals maps at the three frequencies are not fully correlated. Further work (e.g., Vacher et al. 2022a) is needed to model theoretically the impact of polarization angle variations on *EE* and *BB* power spectra of residual maps, which we expect to depend on the correlation between the dust emission properties and the structure of the magnetized interstellar medium.

Our analysis of *Planck* data brings out significant differences between the dust polarization *EE* and *BB* SEDs, and with respect to total intensity, setting new requirements for simulations of the dust polarized foreground and component separation methods. A significant refinement to dust modeling is necessary to ensure an unbiased detection of CMB primordial *B*-modes at the precision required by future CMB experiments.

**Acknowledgements.** A.R. acknowledges financial support from the French space agency (Centre National d'Etudes Spatiales, CNES) and the Italian Ministry of University and Research - Project Proposal CIR01_00010. F.B. acknowledges support from the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022) and CNES. The authors would like to thank the anonymous referee for a thorough reading of the article and for all suggestions and comments that significantly improved the understanding of the text.

## References

Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv e-prints [arXiv:1610.02743]

Abitbol, M. H., Hill, J. C., & Johnson, B. R. 2016, MNRAS, 457, 1796

Ade, P. A. R., Aguirre, J., Ahmed, Z., et al. 2019, JCAP, 2019, 056

Ade, P. A. R., Ahmed, Z., Amiri, M., et al. 2022, ApJ, 927, 77

Ashton, P. C., Ade, P. A. R., Angilè, F. E., et al. 2018, ApJ, 857, 10

Azzoni, S., Abitbol, M. H., Alonso, D., et al. 2021, JCAP, 2021, 047

Benabed, K., Cardoso, J. F., Prunet, S., & Hivon, E. 2009, MNRAS, 400, 219

Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20

Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195

Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914

de la Hoz, E. 2022, ArXiv e-prints [arXiv:2203.04861]

Delouis, J. M., Pagano, L., Mottet, S., Puget, J. L., & Vibert, L. 2019, A&A, 629, A38

Désert, F.-X. 2022, A&A, 659, A70

Fixsen, D. J., Bennett, C. L., & Mather, J. C. 1999, ApJ, 526, 207

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16

Guth, A. H. 1981, Phys. Rev. D, 23, 347

Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1

Hensley, B. S., & Bull, P. 2018, ApJ, 853, 127

Hensley, B. S., & Draine, B. T. 2022, ApJ, submitted [arXiv:2208.12365]

Ichiki, K., Kanai, H., Katayama, N., & Komatsu, E. 2019, Prog. Theor. Exp. Phys., 2019, 033E01

Kamionkowski, M., & Kovetz, E. D. 2016, ARA&A, 54, 227

Krachmalnicoff, N., Carretti, E., Baccigalupi, C., et al. 2018, A&A, 618, A166

Linde, A. 1982, Phys. Lett. B, 116, 335

LiteBIRD Collaboration (Allys, E., et al.) 2022, ArXiv e-prints [arXiv:2202.02773]

Mangilli, A., Aumont, J., Rotti, A., et al. 2021, A&A, 647, A52

Martire, F. A., Barreiro, R. B., & Martínez-González, E. 2022, JCAP, 2022, 003

Mennella, A., Bersanelli, M., Butler, R. C., et al. 2011, A&A, 536, A3

Osumi, K., Weiland, J. L., Addison, G. E., & Bennett, C. L. 2021, ApJ, 921, 175

Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16

Planck Collaboration IX. 2014, A&A, 571, A9

Planck Collaboration XI. 2014, A&A, 571, A11

Planck Collaboration Int. XVII. 2014, A&A, 566, A55

Planck Collaboration Int. XXII. 2015, A&A, 576, A107

Planck Collaboration VIII. 2016, A&A, 594, A8

Planck Collaboration X. 2016, A&A, 594, A10

Planck Collaboration Int. XXX. 2016, A&A, 586, A133

Planck Collaboration Int. XLVI. 2016, A&A, 596, A107

Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109

Planck Collaboration Int. L. 2017, A&A, 599, A51

Planck Collaboration III. 2020, A&A, 641, A3

Planck Collaboration IV. 2020, A&A, 641, A4

Planck Collaboration V. 2020, A&A, 641, A5

Planck Collaboration VI. 2020, A&A, 641, A6

Planck Collaboration XI. 2020, A&A, 641, A11

Planck HFI Core Team 2011, A&A, 536, A4

Puglisi, G., Fabbian, G., & Baccigalupi, C. 2017, MNRAS, 469, 2982

Rosset, C., Tristram, M., Ponthieu, N., et al. 2010, A&A, 520, A13

Svalheim, T. L., Andersen, K. J., Aurlien, R., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202243160

Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90

Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821

Vacher, L., Aumont, J., Boulanger, F., et al. 2022a, A&A, submitted [arXiv:2210.14768]

Vacher, L., Aumont, J., Montier, L., et al. 2022b, A&A, 660, A111

Vacher, L., Chluba, J., Aumont, J., Rotti, A., & Montier, L. 2023, A&A, 669, A5

Weiland, J. L., Addison, G. E., Bennett, C. L., Halpern, M., & Hinshaw, G. 2022, ApJ, 936, 24

Zonca, A., Thorne, B., Krachmalnicoff, N., & Borrill, J. 2021, J. Open Source Softw., 6, 3783