© A. Mangilli et al. 2021

1. Introduction

The precise characterization of the properties of the polarized dust emission from our Galaxy is a crucial issue in the quest for the primordial B-modes of the cosmic microwave background (CMB). If measured, this faint cosmological signal imprinted by the primordial gravitational wave background would be evidence of the inflation epoch and could be used to quantify its energy scale, providing a rigorous test of fundamental physics far beyond the reach of accelerators (Polnarev 1985; Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997). However, accurate determination of diffuse CMB B-mode foregrounds – among which the polarized Galactic dust emission dominates at observing frequencies ≳70 GHz (see e.g., Krachmalnicoff et al. 2016; Planck Collaboration XI 2020) – is required to get an unbiased estimate of the ratio r between tensor and scalar primordial perturbations, a parameter of unknown amplitude scaling the CMB B-mode power on the sky and directly linked to the energy scale at which inflation occurred.

The frequency dependence of the dust emission, assessed through its spectral energy distribution (SED) is one of the characteristics that needs to be determined with the highest accuracy. The Planck data have shown that the SED of the dust emission for total intensity and polarization can be fitted by a modified black body (MBB) law (Planck Collaboration XXX 2014; Planck Collaboration XI 2014; Planck Collaboration Int. XXII 2015; Planck Collaboration XI 2020). Maps of the dust MBB spectral indices and temperatures have been fitted to the total intensity Planck data (see e.g., Planck Collaboration Int. XVII 2014; Planck Collaboration X 2016). These provide evidence that dust emission properties vary across the sky. The data do not provide comparable observational evidence for polarization due to insufficient signal to noise ratio (S/N). However, analyzing total intensity data is directly relevant to polarization. Indeed, Planck and the Balloon-borne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) observations have shown that the polarization fraction is remarkably constant from the far-infrared to millimeter wavelengths, suggesting that the polarized and total dust emission arise predominantly from a single grain population (Planck Collaboration Int. XXII 2015; Gandilo et al. 2016; Ashton et al. 2018; Guillet et al. 2018; Shariff et al. 2019; Planck Collaboration XI 2020).

In the inference of cosmological parameters from CMB data, the spectral frequency dependence of dust polarization and its angular structure on the sky are most often assumed to be separable. This is a simplifying assumption that needs to be improved upon. When spatial variations of the dust emission law are present but ignored, biases that compromise the cosmological analysis can arise and bias the sought-after CMB B-mode signal (Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017; Poh & Dodelson 2017). A bias can also be introduced by additional Galactic emission components or if the dust emission is not perfectly fitted by a MBB emission law. Even if the MBB is observed to be a good fit to the data, dust models do anticipate such departures (e.g., Draine & Hensley 2013; Hensley & Bull 2018).

Spatial SED variations induce a decorrelation between dust maps in different frequency bands, causing a loss of power in the map cross-correlation compared to the geometrical mean of power spectra, which can generate a bias in CMB spectra if it is not properly accounted for. A tentative detection of this effect with the Planck polarization dust data (Planck Collaboration Int. L 2017) was later dismissed (Sheehy & Slosar 2018). This analysis was further extended by Planck Collaboration XI (2020), performing a global multi-frequency fit of the polarized Planck HFI channels with a spectral model that includes decorrelation. These latter authors derive an upper limit on frequency decorrelation that depends on an ad-hoc, possibly misleading model. This shortcoming also applies to the analysis of the latest BICEP2/Keck data (Keck Array & BICEP2 Collaborations 2018). Although current data analyses do not provide conclusive evidence for frequency decorrelation, this is an effect that must be present at some level. The spectral modeling of polarized dust emission has thus become one of the main challenges in the quest of primordial B-modes.

An additional complexity is due to the fact that in sky maps, local variations in dust emission properties across the interstellar medium are averaged along the line of sight and within the beam. When computing power spectra, further averaging occurs through the computation of the expansion of spherical harmonics. Even if the emission law of the dust in any infinitesimal volume element of the Galaxy is perfectly described by a MBB law, after averaging, the dust spectral frequency dependence is no longer a MBB and SED distortions arise. The dust frequency dependence and its angular structure on the sky become interdependent in ways that are difficult to model and generally involve nonlinear transformations.

Chluba et al. (2017) proposed a way to describe the variations of the spectral properties along the line of sight, inside the beam, and across the sky using the moment decomposition around the MBB in the map pixel space. The moments can capture SED distortions due to variations in the dust temperature and spectral index, along the line of sight or between lines of sight, and also in the potential contribution of minor dust emission components. The present paper extends the moment formalism from the map to the angular cross-power spectra, which are highly relevant to CMB B-mode analyses, and to getting rid of noise bias and uncorrelated systematic effects. Similar extensions could also play an important role in the extraction of primordial CMB spectral distortions that are caused by energy release (Zeldovich & Sunyaev 1969; Sunyaev & Zeldovich 1970b,a; Chluba & Sunyaev 2012) and which could be targeted in the future (Kogut et al. 2019; Chluba et al. 2019).

We present the formalism and assess its ability to fit simulations of increasing complexity before using it to analyze Planck High Frequency Instrument (Planck-HFI) total intensity data and build a spectral model in terms of harmonic space moments of the dust spectral index. While the polarization data available today are not sensitive enough to perform such an analysis, Planck-HFI total intensity data offer the required sensitivity and frequency coverage to build a direct spectral model based on astrophysical grounds. Here, we consider Planck intensity data as a proxy to data from future CMB B-modes experiments from Space (Hazumi et al. 2019; Hanany et al. 2019), in terms of S/N and frequency coverage.

The paper is organized as follows. In Sect. 2, we describe the formalism of the dust moment expansion in harmonic space including angular cross-power spectra. Section 3 details our methodology and implementation of the dust moments analysis. We present the simulations and the actual Planck data to which we fit our spectral model in Sect. 4; the results of the fits are presented in Sect. 5. In Sect. 6, within a simplified framework, we relate our spectral model to the analysis of CMB B-modes data. We summarize the main results and present our conclusions in Sect. 7. The paper has five Appendices. The first three detail the simulations (Appendix A), the cross-spectra covariance matrix (Appendix B), and additional fit results (Appendix C). In Appendix D we asses the impact of synchrotron emission on our analysis and in Appendix E we assess the effect of the spatial variations of dust temperatures.

2. Formalism

In this section we present the formalism to describe the moment expansion of the dust intensity SED built from angular cross-power spectra of spherical sky maps. As presented in Chluba et al. (2017), this formalism is a powerful tool with which to account for SED distortions arising from the various averaging effects. We first reiterate the usual and the moment-expansion dust SED parameterizations in Sect. 2.1. In Sects. 2.2 and 2.3 we present the spherical harmonics and the cross-power spectra moment expansion that is used in the analysis we present.

We present the formalism to describe the spectral departures from the MBB associated with derivatives with respect to the dust spectral index β, which is known to vary across the sky (e.g., Planck Collaboration X 2016). This can be easily generalized to include derivatives with respect to the dust temperature (Chluba et al. 2017). However, in the following analysis, we only use the spectral index moment expansion in the Rayleigh-Jeans regime, as temperature and β-variations can have similar effects on the dust SED built from microwave and submillimeter data.

2.1. Dust SED parametrization in the map domain

2.1.1. Modified black body formalism

The commonly used parametrization for a single-temperature dust spectrum is a MBB emission law. We first consider the MBB parametrization without any spatial variations of the SED, that is, with a constant spectral index β₀ and a constant temperature T₀ across the sky. At a given frequency ν the dust intensity map $I_{\text{D}}(\nu ,\widehat{\mathit{n}})$ takes the form:

$$\begin{array}{c}\hfill I_{\mathrm{D}}(\nu ,\widehat{\mathit{n}})=\frac{I_{\nu }(T_0,{\beta }_0)}{I_{{\nu }_0}(T_0,{\beta }_0)}{A}_{\mathrm{D}}(\widehat{\mathit{n}})={\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_0}\frac{B_{\nu }(T_0)}{B_{{\nu }_0}(T_0)}{A}_{\mathrm{D}}(\widehat{\mathit{n}}),\end{array}$$(1)

where $A_{\text{D}}(\widehat{\mathit{n}})$ is in this case a frequency-independent dust intensity template map, ν₀ a reference frequency, and B_ν(T₀) is the Planck law for the temperature T₀. As long as the MBB with constant temperature and spectral index is the correct emission law for all lines of sight, the spatial and the frequency dependence are separable. Nevertheless, different lines of sight probe Galactic regions with very different physical conditions (temperature, dust composition) in the three dimensions of space (we note here that spatial variations between lines of sight and 3D variations along one line of sight lead to similar effects on the final SED). Emission from these regions are averaged along the line of sight and within instrumental beams so that even if the SED of every infinitesimal volume element were to accurately describe by Eq. (1), the MBB would no longer accurately describe the effective emission law. In the following, we formally consider only the variations from one line of sight to another, so that the effective spectral index β or temperature T in Eq. (1) can vary spatially.

2.1.2. Modified black body with spatially varying spectral index

One general attempt to describe the spatial variations of the dust SED in the Rayleigh-Jeans regime of the dust SED is to allow the MBB spectral index to vary spatially. As a consequence, the frequency and spatial dependence of the dust emission are no longer trivially separable. Therefore, the standard MBB formalism of Sect. 2.1.1 must be extended. In this case, the dust intensity map can be written as follows (assuming that the dust temperature remains constant across the sky):

$$\begin{array}{c}\hfill I_{\mathrm{D}}(\nu ,\widehat{\mathit{n}})=\frac{I_{\nu }(T_0,\beta (\widehat{\mathbf{n}}))}{I_{{\nu }_0}(T_0,\beta (\widehat{\mathbf{n}}))}{A}_{\mathrm{D}}(\widehat{\mathit{n}}),\end{array}$$(2)

where now $\beta (\widehat{\mathit{n}})$ accounts for the fact that the spectral index varies from pixel to pixel over the sky. We note that for the sake of clarity, in writing this equation, we ignore the variations along the line of sight.

2.1.3. Modified black body spectral index moment expansion

A more general and powerful parametrization of the dust SED was proposed by Chluba et al. (2017). It introduces the moment expansion of the dust SED, a perturbative expansion of the SED by means of derivatives of the MBB. In our frequency regime, we use derivatives with respect to the dust spectral index β so that the dust map $I_{\text{D}}(\nu ,\widehat{\mathit{n}})$ at a frequency ν now reads:

$$\begin{array}{cc}& I_{\mathrm{D}}(\nu ,\widehat{\mathit{n}})\simeq \frac{I_{\nu }(T_0,{\beta }_0)}{I_{{\nu }_0}(T_0,{\beta }_0)}[A_{\mathrm{D}}(\widehat{\mathit{n}})+{\omega }_1(\widehat{\mathit{n}})ln\left(\frac{\nu }{{\nu }_0}\right)+\frac{1}{2}{\omega }_2(\widehat{\mathit{n}}){ln}^2\left(\frac{\nu }{{\nu }_0}\right)\hfill \\ \hfill & +\frac{1}{6}{\omega }_3(\widehat{\mathit{n}}){ln}^3\left(\frac{\nu }{{\nu }_0}\right)+\cdots ],\hfill \end{array}$$(3)

where ${\omega }_i(\widehat{\mathit{n}})=A_{\text{D}}(\widehat{\mathit{n}})\mathrm{\Delta }\beta {(\widehat{\mathit{n}})}^i$ is the ith order moment map associated to the i^th derivative of the MBB with respect to β (here in an expansion up to third order¹), and $\mathrm{\Delta }\beta (\widehat{\mathit{n}})=\beta (\widehat{\mathit{n}})-{\beta }_0$.

2.2. Dust SED parametrization in spherical harmonics

By definition, decomposition of the dust map into spherical harmonics coefficients implies averaging over various lines of sight over the sky, which is mathematically equivalent to averaging pixels with different SEDs along the line of sight or in an instrumental beam (Chluba et al. 2017). Therefore, we can use the spectral moment decomposition described in Eq. (3). This leads to the following expansion in the spherical harmonics space:

$$\begin{array}{cc}& {(I_{\mathrm{D}})}_{\ell m}^{\nu }=\frac{I_{\nu }(T_0,{\beta }_0(\ell ))}{I{\nu }_0(T_0,{\beta }_0(\ell ))}\times [{(A_{\mathrm{D}})}_{\ell m}+{({\omega }_1)}_{\ell m}ln\left(\frac{\nu }{{\nu }_0}\right)\hfill \\ \hfill & +\frac{1}{2}{({\omega }_2)}_{\ell m}{ln}^2\left(\frac{\nu }{{\nu }_0}\right)+\frac{1}{6}{({\omega }_3)}_{\ell m}{ln}^3\left(\frac{\nu }{{\nu }_0}\right)+\cdots ],\hfill \end{array}$$(4)

where²β₀(ℓ) refers to the effective value of the dust spectral index β₀ for each multipole, as we see in the following. The variables (ω_i)_ℓm (i ∈ {1, …, N}) are the spherical harmonic coefficient of the ith order moment map ${\omega }_i(\widehat{\mathit{n}})$. We note that the spherical harmonics moments (ω_i)_ℓm are not the same as the spatial moments ${\omega }_i(\widehat{\mathit{n}})$, as they involve different averaging.

We stress that the formalism in Eq. (4) accounts for SED distortions in the most general case and not only the ones associated to the averages due to the spherical harmonics decomposition. Therefore, no further extension is required to capture the effects of line-of-sight or beam averaging. We also stress that the line-of-sight and beam-averaging effects, for real-world experiments, can never be avoided, such that $\beta (\widehat{\mathit{n}})$ should already be interpreted as an averaged dust-amplitude-weighted quantity.

2.3. Dust SED parametrization in the cross-power spectra

Based on Eq. (4), we can compute the cross-spectrum between the maps observed in the frequency bands ν₁ and ν₂. Up to third order in terms of β derivatives, this takes the form:

$$\begin{array}{cc}& {\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_2}=\frac{I_{{\nu }_1}(T_0,{\beta }_0(\ell ))I_{{\nu }_2}(T_0,{\beta }_0(\ell ))}{I_{{\nu }_0}^2(T_0,{\beta }_0(\ell ))}\times \{{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{A}_{\mathrm{D}}}\hfill \\ \hfill & 1.\mathrm{o}rder\{\begin{array}{cc}& +[ln\left(\frac{{\nu }_1}{{\nu }_0}\right)+ln\left(\frac{{\nu }_2}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{\omega }_1}\hfill \\ \hfill & +[ln\left(\frac{{\nu }_1}{{\nu }_0}\right)ln\left(\frac{{\nu }_2}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{{\omega }_1{\omega }_1}\hfill \end{array}\hfill \\ \hfill & 2.\mathrm{o}rder\{\begin{array}{cc}& +\frac{1}{2}[{ln}^2\left(\frac{{\nu }_1}{{\nu }_0}\right)+{ln}^2\left(\frac{{\nu }_2}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{\omega }_2}\hfill \\ \hfill & +\frac{1}{2}[ln\left(\frac{{\nu }_1}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_2}{{\nu }_0}\right)+ln\left(\frac{{\nu }_2}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_1}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{{\omega }_1{\omega }_2}\hfill \\ \hfill & +\frac{1}{4}\left[{ln}^2\left(\frac{{\nu }_1}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_2}{{\nu }_0}\right)\right]{\mathcal{D}}_{\ell }^{{\omega }_2{\omega }_2}\hfill \end{array}\hfill \\ \hfill & 3.\mathrm{o}rder\{\begin{array}{cc}& +\frac{1}{6}[{ln}^3\left(\frac{{\nu }_1}{{\nu }_0}\right)+{ln}^3\left(\frac{{\nu }_2}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{\omega }_3}\hfill \\ \hfill & +\frac{1}{6}[ln\left(\frac{{\nu }_1}{{\nu }_0}\right){ln}^3\left(\frac{{\nu }_2}{{\nu }_0}\right)+ln\left(\frac{{\nu }_2}{{\nu }_0}\right){ln}^3\left(\frac{{\nu }_1}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{{\omega }_1{\omega }_3}\hfill \\ \hfill & +\frac{1}{12}[{ln}^2\left(\frac{{\nu }_1}{{\nu }_0}\right){ln}^3\left(\frac{{\nu }_2}{{\nu }_0}\right)+{ln}^2\left(\frac{{\nu }_2}{{\nu }_0}\right){ln}^3\left(\frac{{\nu }_1}{{\nu }_0}\right)]{\mathcal{D}}_{\ell }^{{\omega }_2{\omega }_3}\hfill \\ \hfill & +\frac{1}{36}\left[{ln}^3\left(\frac{{\nu }_1}{{\nu }_0}\right){ln}^3\left(\frac{{\nu }_2}{{\nu }_0}\right)\right]{\mathcal{D}}_{\ell }^{{\omega }_3{\omega }_3}\hfill \end{array}\hfill \\ \hfill & +\cdots \},\hfill \end{array}$$(5)

where the moment cross-power spectra, ${\mathcal{D}}_{\ell }^{\mathit{ab}}$, {a, b}∈{A_D, ω₁, ω₂, ω₃} are defined as³:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\ell }^{\mathit{ab}}=\frac{\ell (\ell +1)}{2\pi }{\displaystyle \sum _{-\ell \le m,m^{\prime }\le \ell }{(a)}_{\ell m}{(b)}_{\ell m^{\prime }}.}\end{array}$$(6)

In Eq. (5), we group terms according to the maximal derivative order in β that appears. These spectral moment truncations are motivated by first constructing the moment maps and then computing all the corresponding cross power spectra, though in this work we measure these cross moment spectra directly from the cross frequency data power spectra. Using this ordering, truncating after the first moment order means including the first three terms, truncating after the second moment order means including the first six terms, and so on. As it is useful in the following, we define the moment functions for a given cross-spectrum, ${\mathcal{M}}_{\mathcal{l}}^{ab}({\nu }_i,{\nu }_j)$, as the moments ${\mathcal{D}}_{\ell }^{\mathit{ab}}$ normalized to the dust amplitude spectrum ${\mathcal{D}}_{\ell }^{A_{\text{D}}{A}_{\text{D}}}$ and re-scaled by the corresponding frequency-dependent numerical factors for the ν_i × ν_j cross-spectrum, as:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\ell }^{\mathit{ab}}({\nu }_i,{\nu }_j)=c^{\mathit{ab}}({\nu }_i,{\nu }_j,{\nu }_0){\mathcal{D}}_{\ell }^{\mathit{ab}}/{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{A}_{\mathrm{D}}},\end{array}$$(7)

where c^ab(ν_i, ν_j, ν₀) are the numerical coefficients which involve the sum and/or the product of the ln(ν_i/ν₀) and ln(ν_j/ν₀) terms, as appearing in Eq. (5). In this way, the ${\mathcal{M}}_{\mathcal{l}}^{ab}$ functions of Eq. (7) show, for a given cross-spectrum, the effective contribution of each moment to the departure of the standard MBB SED. We stress that, in Eq. (5), the dust spectral index β₀(ℓ) is fixed and assigned with optimized values (these steps are described further) for each multipole ℓ.

The model for the dust SED in Eq. (5) can be truncated at any order, depending on the complexity one wants to capture or the available degrees of freedom in the data to be modeled. At zero order (keeping only the first element of the sum, ${\mathcal{D}}_{\ell }^{A_{\text{D}}{A}_{\text{D}}}$), Eq. (5) becomes the cross-power spectrum with a ℓ-dependent spectral index. The higher order terms account for the dust SED distortions to the MBB, meaning that Eq. (5) provides us with a consistent and robust description of the dust SED distortions with respect to the MBB emission law.

Furthermore, as this formalism describes the corrections to the MBB in the angular power spectrum domain – as a function of the multipole ℓ – it naturally provides an efficient framework with which to characterize the frequency decorrelation of the cross-power spectra due to spatial variations of the SED. In contrast to previous attempts (e.g., Planck Collaboration Int. L 2017; Vansyngel et al. 2017; Keck Array & BICEP2 Collaborations 2018), the frequency decorrelation is not introduced by an ad-hoc choice.

We also note that the normalized first-order term can be interpreted as a correction to the MBB spectral index β₀(ℓ) needed to recover the “true” β(ℓ) when spatial and/or line-of-sight variations are present. We therefore define:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\beta }_0(\ell )\equiv {\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{\omega }_1}/{\mathcal{D}}_{\ell }^{A_{\mathrm{D}}{A}_{\mathrm{D}}},\end{array}$$(8)

where the Δβ₀(ℓ) pivot parameter solution depends on the total number of moments that are included and becomes unbiased in the limit of many moments. In fact, this term quantifies the scale-dependent bias that arises from neglecting SED corrections. Adding the moment expansion allows us to eliminate the bias while having nonzero moments at first and higher orders that include the full covariance between the parameters. This connection can be used to introduce an iterative scheme to the analysis, as discussed in Sect. 3.

We stress that β₀(ℓ) in Eq. (8) is related to a power-spectrum weighted average of $\beta (\widehat{\mathit{n}})$ (which, strictly-speaking, itself is a line-of-sight-averaged quantity). A similar averaging process was discussed for the temperature power spectrum stemming from the relativistic Sunyaev-Zeldovich effect (Remazeilles et al. 2019). The higher order moment terms in Eq. (5) quantify (3D) cross-correlations between the spectral indices (Chluba et al. 2017). In general, they each come with their own spatial dependence and maps of these moments can help identify parts of the sky with substantial SED variations.

In summary, Eq. (5) provides, for the first time, a physically motivated model for dealing with spatial and line-of-sight variations of the dust spectral index β at the power-spectrum level.

3. Methodology and implementation

In this section, we detail the methodology and the implementation of the dust moments analysis. The analysis consists of an ℓ-by-ℓ SED fit of a cross-spectra data vector 𝒟_ℓ, gathering all the cross spectra that can be computed from a set of maps M_νi at frequencies ν_i,  i ∈ {1, …, N}:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\ell }\equiv \left[\begin{array}{c}{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1}\\ {\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_2}\\ {\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_3}\\ ⋮\\ {\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N}\\ ⋮\\ {\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N}\end{array}\right]\equiv \left[\begin{array}{c}{\mathcal{D}}_{\ell }(M_{{\nu }_1}\times M_{{\nu }_1})\\ {\mathcal{D}}_{\ell }(M_{{\nu }_1}\times M_{{\nu }_2})\\ {\mathcal{D}}_{\ell }(M_{{\nu }_1}\times M_{{\nu }_3})\\ ⋮\\ {\mathcal{D}}_{\ell }(M_{{\nu }_1}\times M_{{\nu }_N})\\ ⋮\\ {\mathcal{D}}_{\ell }(M_{{\nu }_N}\times M_{{\nu }_N})\end{array}\right].\end{array}$$(9)

The analysis is divided in three steps:

– Step 1: We fit for each multipole ℓ the two parameters β₀(ℓ) and ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$ (T₀ is fixed). Zero-order or MBB fits, in the following, refer to this first step

– Step 2: We fix β₀(ℓ) and then fit Eq. (5). Three parameters are fitted at the first order, six at the second order and ten at the third order.

– Step 3: We update the value of the spectral index to ${\beta }_0^{\text{corr}}(\mathcal{l})={\beta }_0(\mathcal{l})+\mathrm{\Delta }{\beta }_0(\mathcal{l})$, fix it and redo the fits as in step 2. In practice, step 3 may need to be repeated in order to ensure that Δβ₀(ℓ) and thus ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_1}$ become compatible with zero, indicating the convergence.

For instance, in the case of the dust moment expansion up to the third order, the fitted parameters in step 3, for each multipole bin, are: ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$, ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_1}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$, ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_2{\omega }_2}$, ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_3}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_1{\omega }_3}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_2{\omega }_3}$, ${\mathcal{D}}_{\mathcal{l}}^{{\omega }_3{\omega }_3}$. This “three-step” method for the dust moments fit allows us to consistently search for spectral departures from the standard MBB, quantifying, at each multipole, the corrections to the β₀(ℓ) due to SED averaging effects, and therefore providing a description of the scale dependence of the frequency decorrelation.

In practice, we perform a Levenberg-Marquardt least-squares minimization of r = y − m_i, where, at each multipole bin, y ≡ 𝒟_ℓ is the input cross-spectra data vector of Eq. (9) and ${\mathit{m}}_{\mathrm{i}}\equiv {\mathcal{D}}_{\ell }^D$ is the model vector:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\ell }^D=\left[\begin{array}{c}{\mathcal{D}}_{\ell }^{\mathrm{D},{\nu }_1\times {\nu }_1}\\ {\mathcal{D}}_{\ell }^{\mathrm{D},{\nu }_1\times {\nu }_2}\\ {\mathcal{D}}_{\ell }^{\mathrm{D},{\nu }_1\times {\nu }_3}\\ ⋮\\ {\mathcal{D}}_{\ell }^{\mathrm{D},{\nu }_1\times {\nu }_N}\\ ⋮\\ {\mathcal{D}}_{\ell }^{\mathrm{D},{\nu }_N\times {\nu }_N}\end{array}\right].\end{array}$$(10)

For step 1, the model refers to the dust cross-spectra function truncated at zero order, m_i ≡ m₀, that is, the MBB emission with constant temperature and spectral index β₀(ℓ). For steps 2 and 3, the model refers to the dust moment cross-spectra function of Eq. (5), truncated at orders 1, 2, and 3.

In order to take into account the correlations between the cross-spectra 𝒟^{ν_i × ν_j} and 𝒟^{ν_k × ν_l}, the cross-spectra covariance matrix ℂ is included in the fit:

$$\begin{array}{cc}& \mathbb{C}\equiv {\mathcal{C}}_{\mathit{ijkl}}\left(\ell \right)\equiv \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_i\times {\nu }_j},{\mathcal{D}}_{\ell }^{{\nu }_k\times {\nu }_l})\hfill \\ \hfill & =\left[\begin{array}{ccccc}\mathrm{var}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1})& \cdots & \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N})& \cdots & \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1},{\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N})\\ ⋮& \ddots & ⋮& ⋮\\ \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1})& \cdots & \mathrm{var}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N})& \cdots & \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N})\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_1})& \cdots & \mathrm{cov}({\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_1\times {\nu }_N})& \cdots & \mathrm{var}({\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N},{\mathcal{D}}_{\ell }^{{\nu }_N\times {\nu }_N})\end{array}\right].\hfill \end{array}$$(11)

The computation of this cross-spectra covariance matrix from simulations is described in detail, for our applications, in Appendix B. The reduced χ² to be minimized is then defined as:

$$\begin{array}{c}\hfill {\chi }^2=\frac{{\mathit{r}}^T{\mathbb{C}}^{-1}\mathit{r}}{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}},\end{array}$$(12)

where N_d.o.f. is the number of degrees of freedom.

4. Data and simulation settings

In this section we present the Planck intensity data and the simulations used in this paper. We only use the five highest-frequency Planck-HFI channels (143, 217, 353, 545 and 857 GHz). We discard lower frequencies in order to minimize the impact of emission components other than the thermal dust emission which are significant at frequencies lower than 143 GHz but may be neglected at higher frequencies at high Galactic latitude, such as the synchrotron emission, the anomalous microwave emission (AME), and the free-free emission. The CO emission lines at 115, 230, and 345 GHz are significant in the Planck 217 and 353 GHz channels, but their impact can be strongly reduced by applying a tailored mask, which is what we do as detailed later in this section and in Appendix A. The cosmic infrared background (CIB) emission is a significant emission component to the total intensity in all frequency bands from 143 to 857 GHz, with an amplitude relative to dust that increases towards small angular scales (Planck Collaboration XI 2014).

In the following analysis, we consider five different data sets (full sky intensity maps at 143, 217, 353, 545 and 857 GHz): four types of simulations of the Planck intensity data (labeled SIM1, SIM2, SIM3 and SIM4) and the actual Planck intensity data (labeled PR3). All the maps used in this work use the HEALPix⁴ pixelization (Górski et al. 2005).

4.1. Simulated map data sets

For each simulation type, the frequency channel maps $M_{{\nu }_i}^{\text{SIM}}$ are the sum of a noise component N and a sky template S. The sky template S includes a dust component of increasing complexity from SIM1 to SIM4 which has been implemented in order to explore the impact of the spatial variations of the dust spectral index – and eventually other components, such as the CIB – to assess its potential contribution to the moment expansion analysis of the dust.

4.1.1. Noise component

The noise component N is the same for each simulation type. We use 300 Planck End-to-End noise maps ($N_{{\nu }_i}^{\eta },\eta \in \{1,300\}$, 300 realizations per frequency band ν_i) obtained from the FFP10Planck simulations (which include the contributions of both the noise and the residual systematic errors). These maps are publicly available in the Planck Legacy Archive (PLA⁵).

4.1.2. Dust component

The dust component is built from a dust intensity map template at 353 GHz, defined, in MJy sr⁻¹ units, as:

$$\begin{array}{c}\hfill S_{353}^{\mathrm{D}}(\widehat{\mathit{n}})={10}^{20}·M_{{\tau }_{353}}^{\mathrm{D}}(\widehat{\mathit{n}})·B_{\nu =353\mathrm{GHz}}(T_0=19.6\mathrm{K}),\end{array}$$(13)

where B_{ν = 353 GHz}(T₀ = 19.6 K) is the Planck function at 353 GHz for the temperature T₀ = 19.6 K in W m² sr⁻¹, and $M_{{\tau }_{353}}^{\mathrm{D}}(\widehat{\mathit{n}})$ is the dust optical depth map at 353 GHz. For all the simulations, we use the $M_{{\tau }_{353}}^{\text{D}}$ map derived from an MBB fit of Planck dust total intensity maps, obtained with the GNILC component separation method designed to separate dust from CIB anisotropies (Planck Collaboration Int. XLVIII 2016)⁶. The GNILC analysis produced all-sky maps of Planck thermal dust emission at 353, 545, and 857 GHz with reduced CIB contamination. Reducing the CIB contamination of the thermal dust maps is crucial to building maps of dust optical depth, temperature, and dust spectral index that are accurate at high Galactic latitudes. The dust map template $S_{353}^{\mathrm{D}}(\widehat{\mathit{n}})$ at 353 GHz as defined in Eq. (13) is shown in the top panel of Fig. B.2.

The coefficients used to re-scale the dust template at 353 GHz to a different frequency ν_i are defined as:

$$\begin{array}{c}\hfill {\alpha }_{{\nu }_i}(\widehat{\mathit{n}})=\frac{{\nu }_i^{\beta (\widehat{\mathit{n}})}{B}_{{\nu }_i}(T(\widehat{\mathit{n}}))}{(353{\mathrm{GHz})}^{\beta (\widehat{\mathit{n}})}{B}_{353\mathrm{GHz}}(T(\widehat{\mathit{n}}))},\end{array}$$(14)

meaning that the dust map at the frequency ν_i is:

$$\begin{array}{c}\hfill S_{{\nu }_i}^{\mathrm{D}}(\widehat{\mathit{n}})=S_{353}^{\mathrm{D}}(\widehat{\mathit{n}}){\alpha }_{{\nu }_i}(\widehat{\mathit{n}}).\end{array}$$(15)

We define three types of dust simulation, with different levels of complexity in the frequency scaling of the dust intensity in Eqs. (14) and (15):

1. $S_{{\nu }_i}^{{\text{D}}_1}$: simulations with constant dust temperature ($T(\widehat{\mathit{n}})=T_0=19.6$ K) and spectral index ($\beta (\widehat{\mathit{n}})={\beta }_0=1.59$).

2. $S_{{\nu }_i}^{{\text{D}}_2}$: simulations with Gaussian variations of the dust spectral index and fixed temperature ($T(\widehat{\mathit{n}})=T_0=19.6$ K). The Gaussian spectral index map $\beta (\widehat{\mathit{n}})=\mathcal{N}({\beta }_0,\mathrm{\Delta }{\beta }^2)$ is a single random realization (i.e., the same realization for all the simulations) of the normal distribution with a mean β₀ = 1.59 and a 1-σ dispersion Δβ = 0.1. The Δβ value is chosen to roughly match the observed dispersion of the spectral index variations in the Planck data (Planck Collaboration Int. XXII 2015; Planck Collaboration Int. XLVIII 2016). The corresponding $\beta (\widehat{\mathit{n}})$ map is shown in the top panel of Fig. B.3.

3. $S_{{\nu }_i}^{{\text{D}}_3}$: simulations using Planck sky maps of the dust spectral index and temperature. The dust spectral index map $\beta (\widehat{\mathit{n}})$ and the temperature map $T(\widehat{\mathit{n}})$ are those derived from the MBB fit of the Planck GNILC dust total intensity maps (Planck Collaboration Int. XLVIII 2016). For illustration, the $\beta (\widehat{\mathit{n}})$ map is shown in the bottom panel of Fig. B.3.

4. $S_{{\nu }_i}^{{\text{D}}_4}$: simulations using the GNILC Planck sky map of the dust spectral index $\beta (\widehat{\mathit{n}})$ and a constant dust temperature T₀ = 19.6 K. These simulations are not used in the main analysis but presented in Appendix E to assess the effect of the dust temperature spatial variations.

The simulations are computed at the Planck-HFI reference frequencies. For comparison, the PR3 data will be color-corrected to account for the Planck-HFI bandpasses (Planck Collaboration IX 2014).

4.1.3. Cosmic infrared background and synchrotron

In one of our simulation sets, detailed below, we include a CIB component $S_{{\nu }_i}^{\text{CIB}}$. For this, we use the multi-frequency CIB simulation of the Planck Sky Model (PSM) version 1.9 described in Delabrouille et al. (2013)⁷. This CIB map is a random Gaussian realization matching the Planck measured CIB power spectra of Planck Collaboration XVIII (2011). As an illustration, we show the CIB map at 353 GHz in the middle panel of Fig. B.2. Given that the synchrotron emission has a negligible impact on our analysis we do not detail the simulations of this component here, but comment on it in Appendix D.

4.1.4. Simulation types

From the components above, we produce four batches of 100 simulated intensity maps at 143, 217, 353, 545, and 857 GHz, numbered by the superscript η ∈ {1, 100}:

• (a)

  SIM1: $M_{{\nu }_i}^{\text{SIM1},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_1}$,

• (b)

  SIM2: $M_{{\nu }_i}^{\text{SIM2},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_2}$,

• (c)

  SIM3: $M_{{\nu }_i}^{\text{SIM3},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_3}$,

• (d)

  SIM4: $M_{{\nu }_i}^{\text{SIM4},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_3}+S_{{\nu }_i}^{\text{CIB}}$.

We note that for a given simulation type and in a given frequency band, only the noise realization changes from one simulation to another while the sky component remains constant.

4.2. Planck map data set

For the actual Planck map data set, we use the publicly available total intensity maps (observed at 143, 217, 353, 545 and 857 GHz) from the third and latest Planck release. When referring to the data we therefore use the “PR3” label. The CMB component is subtracted from the PR3 intensity data set at the map level. We subtract the SMICA CMB map (Planck Collaboration IV 2020) from the five Planck-HFI frequency channel maps we use in this work. As the SMICA map resolution is different from that of the Planck-HFI frequency channel maps, we first deconvolve the SMICA CMB map from its beam function and then convolve it with the beam of individual Planck-HFI frequency channel maps before subtraction. All the required information to do so is available in the PLA⁵.

4.3. Cross-power spectra computation

The cross-angular power spectra are computed from the data sets (SIM1, SIM2, SIM3, SIM4 and PR3) applying the LR42 sky mask defined in Planck Collaboration Int. XXX (2016) and used in several Planck publications. This mask leaves 42% of the sky for the analysis; it is apodized and includes a galactic mask, a point source mask, and a CO mask, as described in Planck Collaboration Int. XXX (2016), and is shown in Fig. B.1.

To compute the cross-spectra we use the Xpol code described in Tristram et al. (2005). The Xpol code is a pseudo-C_ℓ power spectrum estimator that corrects for the incomplete sky coverage, the filtering effects, and the pixel and beam window functions.

For both the Planck and the simulation data sets, we compute the 15 possible cross-spectra⁸ between the five Planck-HFI channels from 143 to 857 GHz, as described in Sect. 3, and store them in the data vector 𝒟_ℓ (see Eq. (9)). The cross-spectra vector is binned in 15 bins of multipoles⁹ of size Δℓ = 20 in the range most relevant for CMB primordial B-modes analysis (ℓ ∈ {20, 300}), such that 𝒟_b ≡ ⟨𝒟_ℓ⟩_ℓ∈b. As we work only with the binned version of the cross-spectra, in the remainder of this paper and for the sake of clarity, we do not make the distinction between b and ℓ, and so 𝒟_ℓ refers to 𝒟_b.

In order to avoid the noise auto-correlation bias and to reduce the level of correlated systematic errors, we compute the cross-spectra from data split maps (the Planck half-mission maps, HM). We explain how we combine half-mission maps and compute the covariance matrix of the cross-spectra in Appendix B.

After the computation of the cross-spectra 𝒟_ℓ, we subtract, from every data set (PR3 and simulations), the averaged cross-spectrum computed from the 300 Planck End-to-End simulations, which include instrumental noise and systematic effects (see Sect. 4.1.1). This allows us to correct for the small bias linked to residual systematic errors in the data and, by construction, in our simulations. Finally, we apply a color-correction to the PR3 data set cross-spectra to get rid of Planck-HFI-specific calibration effects as described in Planck Collaboration IX (2014). The units of the total intensity cross-spectra for all the data sets are [(MJy sr⁻¹)²].

5. Results for simulations and Planck data

Here we present the results of the fits made for the five data sets (SIM1, SIM2, SIM3, SIM4 and PR3) described in Sect. 4. Each data set is fitted independently for each multipole bin ℓ using the MBB moments expansion with the three-step method introduced in Sect. 3. The maximum order of the moments expansion is limited to third order by the number of degrees of freedom of our data sets (see Appendix B). Here, we focus on our results for the MBB and the moment expansion to this maximum order, but additional plots in Appendix C include results for first and second orders.

This section is organized as follows. We report on the goodness of fits and the dust amplitude spectrum ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$ in Sect. 5.1, on the dust spectral index β₀(ℓ) and its leading order correction Δβ₀(ℓ) in Sect. 5.2, and on higher order moments in Sect. 5.3. Section 5.4 summarizes these results. For each set of simulations, we present mean results averaged over 100 realizations, with error bars corresponding to the standard deviation among realizations. For the PR3 data, we estimate error bars propagating uncertainties through the fit. We checked that this approach, when applied to simulations, yields comparable error bars to those computed from the 100 realizations.

5.1. Goodness of fits and dust amplitude spectra

The goodness of the fits is quantified with the reduced χ²(ℓ) plots, one per data set, shown in Fig. 1. In each plot, we show the reduced χ²(ℓ) for four fits: the MBB law and the moments expansion in Eq. (5) truncated to first, second, and third order.

Fig. 1. Reduced χ2 of the moment expansion fits as a function of the multipole ℓ. From the top left to the bottom right, the panels show the χ2(ℓ) results for the SIM1 (green circles), SIM2 (yellow squares), SIM3 (red diamonds), SIM4 (blue stars), and PR3 (black triangles) data sets. The χ2(ℓ) of the fits at zero (solid black), first (dashed yellow), second (dashed-dotted blue), and third order (dotted red) are displayed.

As expected, the SIM1 set is well fitted by the MBB because these simulations do not include variations of the MBB parameters. The fact that the reduced χ² of the fit has a value close to 1 indicates that our correction for Planck residual systematic errors described in Sect. 4.3 is effective. There is some weak evidence of residual systematic errors in the lowest ℓ = 23 bin; indeed, for this bin a third-order fit is needed for the χ²(ℓ) to reach unity. For SIM2, fitting with higher order terms is required to get a reduced χ²(ℓ)∼1. First order gives a fair χ²(ℓ) (except at very low ℓ) and second order is needed to get a χ²(ℓ) close to unity.

For the SIM3 and SIM4 simulations, as well as for the PR3 data, the MBB yields very poor fits. This illustrates the dual impact of multipole averaging, which introduces spectral complexity and increases the signal-to-noise ratio. For SIM3, a moment expansion up to second-order is required to get a fair fit and up to third order to reach χ²(ℓ)∼1, while for the SIM4 set that includes the CIB, the third order is needed to have a fair reduced χ²(ℓ). For the PR3 data, the reduced χ²(ℓ) is somewhat worse than those for SIM4 (except for the MBB fit which is slightly better), indicating that the Planck data have more spectral complexity than the SIM4 simulations.

Figure 2 shows the amplitude of the dust spectra ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$ for the third-order fit of each data set. We note that these amplitudes do not significantly depend on the order of the fit (see Fig. C.1). As expected, the SIM1, SIM2, and SIM3 amplitudes are essentially indiscernible because they are built from the same dust spatial template and they only differ in the modeling of the dust SED. The SIM4 and PR3 spectra are close to each other. Both depart from the dust-only simulations for multipoles ℓ ≳ 100, that is, for angular scales where the contribution from the CIB component is significant. We note that the power measured on the SIM4 simulations is somewhat larger than that measured for the Planck data for multipoles ranging from about 100 to 220. This may be due to the fact that the GNILC maps used as templates in the dust simulations are not fully free of CIB (Chiang & Ménard 2019).

Fig. 2. Amplitude of the dust spectrum ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$ as a function of the multipole ℓ. The symbols are green circles, yellow squares, red diamonds, and blue stars for the simulation sets SIM1 to SIM4, and black triangles for the Planck PR3 data. The error bars are smaller than the symbols.

5.2. Dust spectral index

The dust spectral indices derived from our fitting are presented in Fig. 3. The top plot displays the spectral indices β₀(ℓ) derived from the MBB fit, and the bottom plot the corrected dust spectral index ${\beta }_0^{\text{corr}}(\mathcal{l})\equiv {\beta }_0(\mathcal{l})+\mathrm{\Delta }{\beta }_0(\mathcal{l})$ in Eq. (8), obtained by fitting the spectral model in Eq. (5) up to third order. The correction depends on the truncation order of the fit as illustrated in Fig. C.2.

Fig. 3. Upper panel: spectral index β0(ℓ) as a function of the multipole ℓ for the MBB fit (step 1). Bottom panel: corrected dust spectral index: ${\beta }_0^{\text{corr}}(\mathcal{l})={\beta }_0(\mathcal{l})+\mathrm{\Delta }{\beta }_0(\mathcal{l})$ from step 3 of the model fit. The symbols for the different data sets are as in Fig. 2.

We find that β₀(ℓ) matches the input spectral index of the simulation β₀ = 1.59 for the SIM1 set. For SIM2, β₀(ℓ) values are slightly larger than the input value at ℓ ≳ 100 but the small difference is corrected when computing ${\beta }_0^{\text{corr}}(\mathcal{l})$.

Our results are less straightforward for the additional sets. For SIM3, we find ${\langle {\beta }_0^{\text{corr}}(\mathcal{l})\rangle }_{\text{SIM3}}=1.537\pm 0.003$ with no systematic dependence with ℓ. For comparison, we computed the spectral index for a MBB from the ratio between the cross-spectra of the 217 and 353 GHz maps and the 353 GHz power spectrum, as in Planck Collaboration XI (2020). We find values in the range 1.52 to 1.55 with no systematic dependence on ℓ, in good agreement with the corrected index from our model. We also computed the median spectral index of the GNILC input map (corrected to the reference temperature T₀ = 19.6 K we use in our model) over the L42 mask. The resulting value, 1.6, computed giving equal weights to each pixel, is slightly larger than ${\langle {\beta }_0^{\text{corr}}(\mathcal{l})\rangle }_{\text{SIM3}}$.

The comparison of SIM3 and SIM4 in Fig. 3 shows that the CIB lowers both β₀(ℓ) and ${\beta }_0^{\text{corr}}(\mathcal{l})$ for ℓ > 100. Above this multipole, the difference between SIM3 and SIM4 values of ${\beta }_0^{\text{corr}}(\mathcal{l})$ increases steadily with ℓ, as does the CIB contribution. The values of ${\beta }_0^{\text{corr}}(\mathcal{l})$ obtained on the PR3 data are systematically lower than the corresponding values for SIM4, but the dependence on ℓ is similar for both sets. We direct the reader to Fig. C.2, which shows a systematic dependence of ${\beta }_0^{\text{corr}}(\mathcal{l})$ on the order of the expansion, that is, on the spectral model. It is interesting to note that from second to third order the values of ${\beta }_0^{\text{corr}}(\mathcal{l})$ move in the opposite direction for SIM4 and PR3. While the two sets of values roughly match for the second-order fits, they differ at third order.

To conclude this section, we compare our PR3 results with earlier determinations of the dust spectral index obtained by analyzing Planck total intensity data.

Planck Collaboration XI (2020) measured the dust spectral index from the ratio between the 217 × 353 and the 353 × 353 cross-spectra over the multipole range ℓ ∈ {4, 170}. For ℓ < 100, where our analysis indicates that the measured spectral index is not biased by the CIB contribution, the spectral indices reported in their Table C.4 for the L42 mask are in the range 1.47–1.50, a little larger than our value ${\langle {\beta }_0^{\text{corr}}(\mathcal{l})\rangle }_{\text{PR3}}=1.45\pm 0.01$ averaged over the same range of multipoles. In Planck Collaboration XI (2020), the spectral index increases with ℓ for ℓ > 100 up to 1.53 ± 0.01 in their ℓ = 150 bin. In our analysis, we observe an opposite trend with the spectral index decreasing for ℓ > 100. It is not clear to us how to interpret this result. We are inclined to consider this as evidence of a spectral mismatch between the SIM4 simulations and the PR3 data, rather than an outcome of the spectral model used to determine the spectral indices. This hypothesis is supported by the fit results in Fig. C.2, which show that the order of the moments expansion does not have a significant effect on the ℓ-dependence.

5.3. High-order moments

In this section, we discuss the results for the high-order (first- to third-order) moments in Eq. (5). We choose to show the moment functions ${\mathcal{M}}_{\mathcal{l}}^{ab}({\nu }_i,{\nu }_j)$ defined in Eq. (7) with ν_i = 143 GHz and ν_j = 545 GHz. We stress that this specific choice of frequencies only impacts the scaling pre-factor c^ab(ν_i, ν_j, ν₀) in Eq. (7), and not the ℓ-dependence of the moments or the relative values between data sets for a given moment. The moment functions are presented in five plots in Fig. 4, one for each data set. All moments in this figure were obtained from fits to third order. The top left panel of Fig. 4 represents the first-order moment function ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_1}$, which is made compatible with zero through iteration on Δβ₀(ℓ) (see Sect. 3).

Fig. 4. Step 3 moment functions ${\mathcal{M}}_{\mathcal{l}}^{ab}$ as a function of the multipole ℓ defined in Eq. (7) for the 143 × 545 cross-spectrum. The plots refer, from top left to bottom right, to ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_1}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$, ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_2}$, ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_3}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_3}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_3}$, and ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_3{\omega }_3}$. The symbols refer to the different data sets: SIM1 (green circles), SIM2 (yellow squares), SIM3 (red diamond), SIM4 (blue stars), and PR3 (black triangles).

As expected, we do not detect any moment function for SIM1. This result gives us confidence that residual systematic errors included in our noise simulations (see Sect. 4.1.1) do not have a significant impact on high-order moments. For the SIM2 data set, the ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$ function is detected with an amplitude increasing with ℓ, while the other moments are consistent with zero. For SIM3, we find that the ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$ moment is nearly scale-independent and is significantly larger than that of SIM2 at low ℓ. For this set, two additional moment functions are detected, namely ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$ and more marginally ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$. Comparing moment functions for the SIM3 and SIM4 sets, we find that the CIB has a significant impact on several moments at ℓ ≳ 100, but not on ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$. Nearly all of the other moment functions increase with ℓ and are close to zero in the lowest ℓ-bins.

For the PR3 data, all the moment functions (from first- to third-order) are detected with an absolute amplitude larger than that measured on the simulated data sets. The ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$ moment function shows a similar amplitude as for SIM3 and SIM4 for ℓ ≳ 100, while at larger angular scales it does not. For the ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$ and ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_2}$, the PR3 data set has an overall ℓ behavior close to that of SIM4 but with an increased absolute value. The ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_3}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_3}$ and ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_3{\omega }_3}$ of PR3 match those of SIM4 for ℓ ≲ 70 but progressively deviate from them for higher multipoles. Finally, we point out that the ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_3}$ moment function is the one with highest absolute amplitude for the PR3 set. Surprisingly, it has the opposite sign, and a different ℓ-dependence from the SIM4 moment.

5.4. Discussion

Here, we summarize the main results of the fits before briefly discussing our interpretation.

– The goodness of the fit obtained for the simulations demonstrates the ability of the moment expansion to account for spatial variations of the dust SED, even when the MBB law provides a very poor fit. Except for the simplest SIM1 simulations, with constant temperature and spectral index, high-order moments are significantly detected for the other simulations and the PR3 data.

– When spatial variations of the dust SED are present, the spectral index inferred from the MBB fit is biased. Fitting high-order moments, we obtain a significantly different value that cancels the first-order moment Aω₁.

– The comparison of the moments obtained for the SIM3 and SIM4 simulations, which only differ by the addition of the CIB, shows that the CIB has a significant impact on moments at ℓ ≳ 100. To account for this additional emission component, we need to extend the moment expansion to third order.

– The moments on the SIM4 simulations that include realistic spatial SED variations and the CIB are quantitatively different from those measured on the PR3 data.

The moment functions in Fig. 4 are difficult to interpret at ℓ > 100 due to the CIB contribution, but for lower multipoles one may relate the results to dust emission properties. For the PR3 data at ℓ < 100, the three most significant moment functions in decreasing order are ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_3}$, ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$ and ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$. We point out that the simulations do not match any of these three moment functions. The most immediate interpretation of this mismatch is the lack of variations of the dust MBB parameters along the line of sight in the simulations, but this may not be the sole explanation. The dust emission could also comprise two or more emission components, which are not fully correlated on the sky (Draine & Hensley 2013; Guillet et al. 2018; Hensley & Bull 2018).

It is not straightforward to provide a specific interpretation of the amplitude, the scale dependence, or the hierarchy of the moments fits. One difficulty lies in the fact that the moments expansion does not decompose the data into independent components: the high-order moments depend on the expansion order as illustrated in Figs. C.3–C.6. Furthermore, the moments also quantify the SED averaging that occurs when going from pixel space to harmonic space. Without a model of the sky emission, it is therefore difficult to link a given high-order moment to a physical property of the dust emission. An iteration on dust and CIB simulations converging towards a moment decomposition in agreement with that of the Planck data would be needed to quantify possible interpretations. We foresee that the moment decomposition could be used as a quantitative metric to obtain simulations that better match the data, but this is beyond the scope of the present work.

6. Implications for the tensor-to-scalar ratio measurement

We finally discuss the potential impact of our results on the measurement of the tensor-to-scalar ratio r. We performed the moment expansion of the dust SED from the Planck intensity power spectra. There is no reason for the departures from the MBB SED that we observed and quantified to be absent in polarization. Here, we indeed assume that the SED of the dust B-modes shows spectral departures from the mean MBB SED of the same relative order as the ones we observed in intensity and we derive the potential bias on r that would result from neglecting them.

In order to be conservative in this determination, we use the SIM3 simulated data set: we know from our χ² analysis of Sect. 5.1 that the actual dust intensity in the Planck PR3 data contains at least the same level of departure from the MBB that is present in this simulated data set. We could have used SIM4 or PR3 data sets, but as we see in Sect. 5, it is not trivial to decipher which part of the SED distortion is due to dust and CIB; the latter component being expected to contribute much less than the dust to the B-modes.

We consider the SIM3 cross-spectra fit with Eq. (5) at different orders in the moment expansion. We look here at the relative difference between the SIM3 data set cross-spectra and the fitted model of Eq. (5), for every cross-spectra between frequencies ν_i and ν_j. We focus on the multipole bin centered at ℓ₀ = 80, as this scale corresponds to the CMB primordial B-modes peak. This relative difference reads:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathcal{D}}_{{\ell }_0}({\nu }_i\times {\nu }_j)\equiv \frac{{\mathcal{D}}_{{\ell }_0}^{\mathrm{SIM}3}({\nu }_i\times {\nu }_j)-{\mathcal{D}}_{{\ell }_0}^{\mathrm{fit}}({\nu }_i\times {\nu }_j)}{{\mathcal{D}}_{{\ell }_0}^{\mathrm{SIM}3}({\nu }_i\times {\nu }_j)}.\end{array}$$(16)

This relative difference is displayed in Fig. 5. Focusing on the 143 × 143 cross-spectrum, which is a frequency channel indicative of typical CMB B-mode experiments, the simple MBB fit leaves a Δ𝒟_ℓ0(143 × 143) = 10.9% residual, the first-order fit leaves Δ𝒟_ℓ0(143 × 143) = 3.2%, the second-order fit leaves Δ𝒟_ℓ0(143 × 143) = 0.5% and the third-order fit leaves Δ𝒟_ℓ0(143 × 143) = 0.06%.

Fig. 5. Upper panel: mean SED in MJy2 sr−1 for the 100 SIM3 cross-spectra as a function of the effective frequency ($\sqrt{{\nu }_i·{\nu }_j}$, in GHz) for the multipole bin centered at ℓ = 80. The MBB (solid black), first- (dashed yellow), second- (dashed-dotted blue), and third-order (dotted red) best fits are displayed and not distinguishable. Bottom panel: relative percentage difference between the SIM3 SED and the MBB (black diamonds) and the first- (yellow diamonds), second- (blue diamonds) and third-order (red diamonds) best fits.

Let us now consider a future CMB B-mode experiment that has the same frequency channels and a S/N for the dust B-modes that is similar to those of Planck-HFI for the dust intensity (e.g., LiteBIRD, Hazumi et al. 2019). We know from Planck Collaboration Int. XXX (2016) that the dust B-mode power spectrum computed on the LR42 mask has an amplitude of A_D(ℓ₀) = 78.6 μK² at 353 GHz. If we convert this amplitude into the r-equivalent amplitude at 150 GHz r_D (Planck Collaboration Int. XXX 2016), it corresponds to r_D = 1.8. Our results therefore suggest that a CMB B-mode experiment looking at half the sky could find an analysis bias of Δr = 0.109 ⋅ r_D = 0.20 by assuming that the dust B-modes follow a MBB SED, Δr = 0.032 ⋅ r_D = 0.06 by assuming an first-order moment expansion, Δr = 0.005 ⋅ r_D = 0.009 by assuming a second-order moment expansion, and Δr = 0.0006 ⋅ r_D = 0.001 by assuming a third-order moment expansion (as we see in the above analysis, parameterizing a dust decorrelation amounts to fitting the first order, if the parametrization is correct).

The region of the sky observed by the BICEP2/Keck experiment has r_D = 0.11 (Keck Array & BICEP2 Collaborations 2018). If we transpose our results to this region, we see that a MBB fit of the dust B-modes would lead to Δr = 0.109 ⋅ r_D = 0.01, a decorrelation or first-order analysis would lead to Δr = 0.032 ⋅ r_D = 4 × 10⁻³, a second-order analysis to Δr = 0.005 ⋅ r_D = 5 × 10⁻⁴, and a third-order analysis to Δr = 0.0006 ⋅ r_D = 7 × 10⁻⁵.

The Δr values in the different cases are summarized in Table 1. Although these Δr values are rough estimates (that might be overestimated in the BICEP2/Keck case because fewer SED spatial variations could occur on this small region), they provide an insight into the order of magnitude of the potential bias. The values seen are in strong support of the need to take into account the spectral departures from the dust MBB in future CMB B-mode analyses, targeting r values down to 10⁻³ and beyond.

These conclusions are further supported by the moment decomposition of the Planck PR3 data set at ℓ < 100, where the CIB contribution is likely to be negligible. As the moment functions ${\mathcal{M}}_{\mathcal{l}}^{ab}$ measure a fractional departure from a simple MBB law, we can see that in order to reach an accuracy in the dust subtraction of 10⁻³, we need to consider all moments with an absolute amplitude larger than 10⁻³. As can be seen in Fig. 4, most of the moments in the PR3 data set decomposition up to third order have an absolute amplitude that is greater than this threshold. In that sense, dust-dominated angular scales of the intensity PR3 data set additionally stress the need for a third-order expansion of the polarized dust SED in order to reach the accuracy targeted by future CMB B-mode experiments.

7. Summary and conclusions

In this paper, we present a model that describes the Galactic dust SED for total intensity at Planck-HFI frequencies in terms of SED distortions with respect to the MBB emission law. This model accounts for variations in the dust SED on the sky and along the line of sight in order to provide an astrophysically motivated description at the power spectrum level. The model formalism relies on expansion of the dust emission SED in moments around the MBB law related to derivatives with respect to the dust spectral index. These high-order moments lead to frequency decorrelation; departures from the MBB inevitably appear because of averaging effects along the line of sight and within the beam, and, as is most relevant here, because of the spherical harmonic expansion performed in the data processing.

We applied our analysis to total intensity cross-spectra computed from the combination of CMB-corrected PR3 Planck data at the five HFI channels at 143, 217, 353, 545, and 857 GHz, and to four sets of foreground simulations of increasing complexity. The main conclusions of our analysis can be summarized as follows.

Our analysis quantifies the spectral complexity of the Planck total intensity data at frequencies larger than 143 GHz. At ℓ ≳ 100, the CIB is a significant component contributing to the complexity. At lower multipoles, the data are dust dominated but the dust simulations based on MBB parameters fitted on Planck maps fail to match the most significant moments. In future work, the moment decomposition could be used to obtain improved simulations of dust and CIB emission, which better match Planck data and may be used to test possible interpretations.

We extend our results to B-mode analyses within a simplified framework. We find that neglecting the dust SED distortions of the dust polarization with respect to the MBB, or trying to model them with a single ad-hoc parameter, could lead to biases larger than the accuracy of the component separation required to search primordial B-modes down to a tensor-to-scalar ratio r = 10⁻³.

If our results extend to polarized emission without any additional complexity, we anticipate that moment expansion up to third order would be required to model the dust polarization SED to the accuracy of future CMB B-mode experiments. If this is a valid statement, it sets constraints on the number of frequency bands required to separate dust and CMB polarization. At least four and five dust-dedicated frequency channels are needed in order to perform second- and third-order moment expansion fits, respectively.

Additional difficulties for B-mode searches could arise from changes in polarization angles across frequencies, which would make the decomposition of polarized dust emission in E and B-modes frequency-dependent. Further complexity may arise because of variation in dust temperature, which we did not include here. Similarly, synchrotron foregrounds at low frequencies will require an independent moment expansion.

Based on these findings, we conclude that the moment expansion of the dust SED is a new promising tool to model the dust component at the level of precision needed for the measurement of the CMB primordial B-modes. This paper presents a first step in this direction, providing the formalism and the first qualitative results based on the Planck total intensity data and simulations. When dealing with the polarization, other details should be carefully considered and added, such as for instance the fact that the magnetic field direction will project variations of the SED differently in Q and U, which implies that, generally, two independent moment expansions are needed. Studying the moment expansion method specifically for polarization will be another important step and the focus of a forthcoming publication.

It will also be important to study applications of the power spectrum moment expansion for extractions of primordial CMB spectral distortions. The expected signals are small (e.g., Chluba 2016) and heavily obscured by foregrounds. Most extraction methods mainly use information based on the SED shapes and neglect spatial information (Sathyanarayana Rao et al. 2015; Desjacques et al. 2015; Abitbol et al. 2017), but a more recent study explores the benefits of using spatial information (Rotti & Chluba 2021). Spatial information could be further exploited using the techniques described here, which warrants further investigation.

Acknowledgments

This research was supported by the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022). A.M. acknowledges the support of the DIA-ASTRO fellowship from the French space agency (Centre National d’Etudes Spatiales, CNES). A.R. is supported by the ERC Consolidator Grant CMBSPEC (No. 725456) as part of the European Union’s Horizon 2020 research and innovation program. J.C. was supported by the Royal Society as a Royal Society University Research Fellow at the University of Manchester, UK.

References

Appendix A: Mask and foreground templates

In this section we provide details on the mask and on the foreground templates that we used to generate the simulations. For both the data and the simulations, we used a mask that combines a 50% apodized galactic cut and a point source mask. This mask, referred to as LR42 and defined in Planck Collaboration Int. XXX (2016), has been extensively used in the Planck analyses, has f_sky = 0.42, and is shown in Fig. B.1.

As described in Sect. 4, we generated different sets of dust and multi-component simulations. We provide an illustrative example of the foreground templates at 353 GHz in MJy sr⁻¹ units in Fig. B.2. From top to bottom, the figure shows the dust template defined in Eq. (13), the CIB template, and the synchrotron template.

Given that investigating the impact of the spatial variations of the dust spectral index is a key part of our analysis, we show in Fig. B.3 the $\beta (\widehat{\mathit{n}})$ maps used in the simulations. The top panel of Fig. B.3 shows the $\beta (\widehat{\mathit{n}})$ map used to generate the SIM2 dust simulations with Gaussian $\beta (\widehat{\mathit{n}})$ variations with Δβ = 0.1 around β₀ = 1.59. The bottom panel of Fig. B.3 shows the GNILC $\beta (\widehat{\mathit{n}})$ map used to generate SIM3 simulations (Planck Collaboration Int. XLVIII 2016).

Appendix B: Cross-spectra definition, covariance, and correlation matrices

Here we provide more details on how we compute the angular power spectra of the data sets we consider in the paper, as introduced in Sect. 4.3, and on how this impacts the cross-power spectra correlations and their statistical independence.

Fig. B.1. LR42 mask used in the analysis.

Fig. B.2. Foreground templates at 353 GHz in MJy sr−1 units. Top panel: dust template defined in Eq. (13). Middle panel: CIB template. Bottom panel: synchrotron template.

Fig. B.3. Dust spectral index map $\beta (\widehat{\mathbf{n}})$ used for the dust simulations SIM2 with Gaussian β variations with Δβ = 0.1 (top panel) and for the dust simulations SIM3 with β variations estimated from the data with the GNILC component separation method (bottom panel).

Fig. B.4. Cross-spectra correlation matrices computed from Eq. (B.2) for the SIM1 data set in the multipole bin centered at ℓ = 100. Upper panel: correlation matrix computed from the cross-spectra as defined in Eq. (B.1). Middle panel: correlation matrix computed for our toy-model of the correlation presented in Appendix B.2. Bottom panel: correlation matrix computed from the cross-spectra as defined in Eq. (B.6).

B.1. Definitions

In order to avoid the noise auto-correlation bias and to reduce the level of correlated systematic errors, we compute the cross-spectra from data split maps (the 2 Planck half-mission maps, HM, namely HM1 and HM2). The general philosophy would be to reproduce Planck Collaboration Int. XXX (2016) to construct the cross-power spectra from two frequency maps M_νi and M_νj:

$$\begin{array}{c}\hfill \{\begin{array}{cc}{{\mathcal{D}}_{\ell }(M_{{\nu }_i}\times M_{{\nu }_j})|}_{i=j}=\hfill & {\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}1}\times M_{{\nu }_i}^{\mathrm{HM}2})\hfill \\ {{\mathcal{D}}_{\ell }(M_{{\nu }_i}\times M_{{\nu }_i})|}_{i\ne j}=\hfill & \frac{1}{4}[{\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}1}\times M_{{\nu }_i}^{\mathrm{HM}1})\hfill \\ \hfill & +{\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}1}\times M_{{\nu }_i}^{\mathrm{HM}2})\hfill \\ \hfill & +{\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}2}\times M_{{\nu }_i}^{\mathrm{HM}1})\hfill \\ \hfill & +{\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}2}\times M_{{\nu }_i}^{\mathrm{HM}2})]\hfill \end{array}\end{array}$$(B.1)

The idea of doing the sum of the four cross-spectra in Eq. (B.1) is to increase the S/N of D_ℓ(M_νi×M_νi)|_i ≠ j, presumably at the cost of the statistical independence between distinct cross-spectra, as we see in the following. To preserve the statistical independence between the cross-spectra, we do not adopt the definition of Eq. (B.1) but that of Eq. (B.6), in Appendix B.3, for the reasons that are presented in Appendix B.2.

The covariance matrix ℂ from Eq. (11), which is used for the fits we perform throughout the paper, is computed from 100 pairs of half-mission maps of a given data set (e.g., SIM1 simulations are used to compute ℂ when dealing with SIM1). As they are the closest to the data in terms of physical components and component complexity, the covariance matrix is inferred from the SIM4 simulations when dealing with the actual Planck data set PR3.

In order to assess the statistical independence between the cross-spectra of our data sets, we build the correlation matrix from the covariance matrix of Eq. (11):

$$\begin{array}{cc}& \mathbb{R}\equiv {\mathcal{R}}_{\mathit{ijkl}}\left(\ell \right)\equiv {\mathrm{corr}}_{\mathit{ijkl}}\left(\ell \right)\equiv \frac{{\mathcal{C}}_{\mathit{ijkl}}\left(\ell \right)}{\sqrt{{\mathcal{C}}_{\mathit{ijij}}\left(\ell \right){\mathcal{C}}_{\mathit{klkl}}\left(\ell \right)}}\hfill \\ \hfill & =\frac{{\mathcal{C}}_{\mathit{ijkl}}\left(\ell \right)}{\sqrt{{\mathrm{var}}_{\mathit{ij}}\left(\ell \right){\mathrm{var}}_{\mathit{kl}}\left(\ell \right)}}=\frac{{\mathcal{C}}_{\mathit{ijkl}}\left(\ell \right)}{{\sigma }_{\mathit{ij}}\left(\ell \right){\sigma }_{\mathit{kl}}\left(\ell \right)}.\hfill \end{array}$$(B.2)

This correlation matrix is displayed for the SIM1 data set in Fig. B.4, for the multipole bin centered at ℓ = 100 (the shape of ℝ is qualitatively the same in each multipole bin). It is significantly nondiagonal, showing large correlations between numerous cross-spectra. This highlights that the cross-spectra, as they are defined in Eq. (B.1), are not all independent. In the following section we see why this happens. Below, we change the definition of the cross-spectra in Eq. (B.1) to minimize these correlations.

B.2. Expected correlations from a toy model

In a high-S/N mixture of interstellar dust and instrumental noise, the correlation between frequency channels cross-spectra can become significant and needs to be taken into account in minimization processes.

Let us suppose that a map of the sky M_i ≡ M_νi observed at a frequency channel ν_i can be written as the sum of a dust component D and an instrumental noise term N so that M_i ≃ D_i + N_i. The cross angular power spectra then read:

$$\begin{array}{c}\hfill M_i\times M_j=D_i\times D_j+D_i\times N_j+D_j\times N_i+N_i\times N_j.\end{array}$$(B.3)

As the dust templates are the same in every simulation, the D × D term does not contribute to the covariance (nor to the variance) and in the high-S/N regime the N × N term can be neglected with respect to the D × N terms. Thus, the correlation coefficient between the power spectra computed from these maps, as defined in Eq. (B.2), reads¹⁰:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathit{ijkl}}\simeq \frac{\mathrm{cov}(D_i\times N_j+D_j\times N_i,D_k\times N_l+D_l\times N_k)}{\sqrt{\mathrm{var}(D_i\times N_j+D_j\times N_i)\mathrm{var}(D_k\times N_l+D_l\times N_k)}}.\end{array}$$(B.4)

Let us now consider as an example the specific correlation between 2 Planck cross-spectra, namely 143 × 217 and 143 × 353 (ν_i = 143, ν_j = 217, ν_k = 143 and ν_l = 353 GHz). As cov(X + Y, Z) = cov(X, Z)+cov(Y, Z), the most significant among the developed terms in the numerator of Eq. (B.4) is cov(D₂₁₇ × N₁₄₃, D₃₅₃ × N₁₄₃) as it is the only one to involve twice the same noise map and hence would be the most “covariant”. If we make the assumption that the dust component is spatially the same at each frequency (i.e., the same dust spatial template D), scaling as D_i = A_i ⋅ D and that the noise basically scales as N_i = B_i ⋅ N (where A_i and B_i are scalars), we find:

$$\begin{array}{cc}& {\mathcal{R}}_{143\times 217,143\times 353}\simeq \hfill \\ \hfill & \frac{A_{217}{A}_{353}{B}_{143}^2}{\sqrt{(A_{217}^2{B}_{143}^2+A_{143}^2{B}_{217}^2)·(A_{353}^2{B}_{143}^2+A_{143}^2{B}_{353}^2)}}.\hfill \end{array}$$(B.5)

For a typical dust MBB SED with β₀ = 1.59 and T₀ = 19.6 K and the Planck-HFI noise levels, we find that ℛ_{143 × 217, 143 × 353} = 0.90, which is a significant correlation. We note that it would be markedly decreased if the noise was independent between the 2 143 GHz maps used in the two different cross-spectra. This can happen using two different data splits (as the Planck HM maps) as we see in the following.

By applying a reasoning equivalent to that of the example above to all the cross-spectra covariances in the correlation matrix ℝ, we can build the full toy-model correlation matrix for the data sets we use in this paper. It is displayed in the middle panel of Fig. B.4. We can see that despite the simplicity of our assumptions, this toy-model correlation matrix qualitatively reproduces that of our data sets (top panel of the same figure).

B.3. Minimizing the correlations and effective degrees of freedom

In order to minimize the correlations between the cross-spectra of our data sets, we change the definition of Eq. (B.1) in the case i ≠ j:

$$\begin{array}{c}\hfill \{\begin{array}{c}{{\mathcal{D}}_{\ell }(M_{{\nu }_i}\times M_{{\nu }_j})|}_{i=j}={\mathcal{D}}_{\ell }(M_{{\nu }_i}^{\mathrm{HM}1}\times M_{{\nu }_i}^{\mathrm{HM}2})\hfill \\ {\mathcal{D}}_{\ell }(M_{143}\times M_{217})\equiv {\mathcal{D}}_{\ell }(M_{143}^{\mathrm{HM}2}\times M_{217}^{\mathrm{HM}1})\hfill \\ {\mathcal{D}}_{\ell }(M_{143}\times M_{353})\equiv {\mathcal{D}}_{\ell }(M_{143}^{\mathrm{HM}1}\times M_{353}^{\mathrm{HM}2})\hfill \\ {\mathcal{D}}_{\ell }(M_{143}\times M_{545})\equiv {\mathcal{D}}_{\ell }(M_{143}^{\mathrm{HM}1}\times M_{545}^{\mathrm{HM}1})\hfill \\ {\mathcal{D}}_{\ell }(M_{143}\times M_{857})\equiv {\mathcal{D}}_{\ell }(M_{143}^{\mathrm{HM}2}\times M_{857}^{\mathrm{HM}2})\hfill \\ {\mathcal{D}}_{\ell }(M_{217}\times M_{353})\equiv {\mathcal{D}}_{\ell }(M_{217}^{\mathrm{HM}2}\times M_{353}^{\mathrm{HM}1})\hfill \\ {\mathcal{D}}_{\ell }(M_{217}\times M_{545})\equiv {\mathcal{D}}_{\ell }(M_{217}^{\mathrm{HM}1}\times M_{545}^{\mathrm{HM}1})\hfill \\ {\mathcal{D}}_{\ell }(M_{217}\times M_{857})\equiv {\mathcal{D}}_{\ell }(M_{217}^{\mathrm{HM}2}\times M_{857}^{\mathrm{HM}2})\hfill \\ {\mathcal{D}}_{\ell }(M_{353}\times M_{545})\equiv {\mathcal{D}}_{\ell }(M_{353}^{\mathrm{HM}1}\times M_{545}^{\mathrm{HM}1})\hfill \\ {\mathcal{D}}_{\ell }(M_{353}\times M_{857})\equiv {\mathcal{D}}_{\ell }(M_{353}^{\mathrm{HM}2}\times M_{857}^{\mathrm{HM}2})\hfill \\ {\mathcal{D}}_{\ell }(M_{545}\times M_{857})\equiv {\mathcal{D}}_{\ell }(M_{545}^{\mathrm{HM}2}\times M_{857}^{\mathrm{HM}1}).\hfill \end{array}\end{array}$$(B.6)

This is the definition of the data sets cross-spectra we adopt in this paper. The correlation matrix built for the SIM1 data set from this definition of the cross-spectra is displayed, for the multipole bin centered at ℓ = 100, in the bottom panel of Fig. B.4. We can see that the correlations have been significantly decreased with respect to those of the cross-spectra defined in Eq. (B.1) and that the correlation matrix is closer from being diagonal.

An eigenvalue analysis of the correlation matrix as computed in Eq. (B.1) indicates that the effective degrees of freedom from the 15 cross-spectra is N_d.o.f. = 5. When defining the cross-spectra as in Eq. (B.6), it becomes N_d.o.f. = 11, allowing us to fit the SED moment expansion of the cross-spectra up to third order (see Sect. 5).

Appendix C: Truncating the dust moment expansion at difference orders

In Sect. 5 we present results from the third-order dust SED moment expansion. Here, we present the corrected dust amplitude ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$, the corrected spectral index ${\beta }_0^{\text{corr}}(\mathcal{l})$, and the relevant moment functions ℳ_ℓ in the case where Eq. (5) is truncated at first and second order in our three-step fitting procedure.

These other truncating order results are displayed in Figs. C.1–C.6 for ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$, ${\beta }_0^{\text{corr}}(\mathcal{l})$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$, ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$, ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$ and ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_2}$, respectively.

Fig. C.1. Corrected dust amplitude ${\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}}$ after step 3 of the fit, when truncating Eq. (5) at first order (“1” label on the plot markers), second order (“2” label), and third order (“3” label, same values as those in Sect. 5). The symbols refer to the different data sets: SIM1 (green circles), SIM2 (yellow squares), SIM3 (red diamond), SIM4 (blue stars), and PR3 (black triangles).

Fig. C.2. Same as Fig. C.1 but for ${\beta }_0^{\text{corr}}(\mathcal{l})$.

Fig. C.3. Same as Fig. C.1 but for ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_1}$, normalized for the 143 × 545 GHz cross-spectrum.

Fig. C.4. Same as Fig. C.1 but for ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$, normalized for the 143 × 545 GHz cross-spectrum.

Fig. C.5. Same as Fig. C.1 but for ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_1{\omega }_2}$, normalized for the 143 × 545 GHz cross-spectrum.

Fig. C.6. Same as Fig. C.1 but for ${\mathcal{M}}_{\mathcal{l}}^{{\omega }_2{\omega }_2}$, normalized for the 143 × 545 GHz cross-spectrum.

SIM1, SIM2, and SIM3 data set results are stable with the truncating order, while SIM4 and PR3 change significantly. Nevertheless, SIM4 and PR3 results are affected in a very different way. For example, ${\beta }_0^{\text{corr}}(\mathcal{l})$ values increase with the truncating order for SIM4, while they evolve in the opposite way for SIM3. However, this behavior is not observed for the moment functions.

Appendix D: Impact of the synchrotron emission

In this section we quantify the impact of the synchrotron emission on the moment analysis by comparing the results between two types of simulations: the SIM4 simulations, described in Sect. 4, and the SIM5 simulations that in addition to SIM4 also include a synchrotron component:

SIM5: $M_{{\nu }_i}^{\text{SIM5},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_3}+S_{{\nu }_i}^{\text{CIB}}+S_{{\nu }_i}^{\text{Sync}}$,

where $S_{{\nu }_i}^{\text{Sync}}$ is the synchrotron template at the frequency ν_i. This component is generated, as for the CIB component, using the Planck Sky Model (PSM) version 1.9 described in Delabrouille et al. (2013). Examples of the CIB and synchrotron templates at 353 GHz are shown in the middle and bottom panels of Fig. B.2, respectively. We can already note that, at this frequency, the synchrotron is sub-dominant with respect to the CIB by at least three orders of magnitude.

Figure D.1 shows the reduced χ² results of the MBB and the dust moment fits up to third order for the SIM4 and SIM5 simulations. The lines are barely distinguishable. The relative difference of the dust amplitude spectrum ${\mathcal{D}}_{\mathcal{l}}^{{\text{A}}_{\text{D}}{\text{A}}_{\text{D}}}$, the dust spectral index β₀(ℓ), the corrected spectral index ${\beta }_0^{\text{corr}}(\mathcal{l})$, and the dust moments functions ${\mathcal{M}}_{\mathcal{l}}^{ab}$ up to third order between SIM4 and SIM5 are shown in Fig. D.2. For most of these quantities, this relative difference is very small. It is bigger for some of them, for example for the ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_i}$ moment functions, but still well within the propagated error bars (due to division by small numbers).

Fig. D.1. Reduced χ2 of the fit of Eq. (5) at zero (solid black), first (dashed yellow), second (dashed-dotted blue), and third order (dotted red) for the SIM5 (including synchrotron, pink plus signs) and SIM4 (blue stars) data sets.

Fig. D.2. Relative difference between SIM5 and SIM4 fitted parameters. ΔX = (XSIM5 − XSIM4)/XSIM4 with $X\in \{{\mathcal{D}}_{\mathcal{l}}^{A_{\text{D}}{A}_{\text{D}}},{\beta }_0,{\beta }_0^{\text{corr}},{\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_1},\text{etc}.\}$ as a function of the multipole ℓ.

According to these results we can therefore conclude that the synchrotron emission has a negligible impact on the dust moment analysis for the Planck-HFI channels from 143 to 857 GHz.

Appendix E: Dust simulations with varying spectral index and constant temperature

In the present work, the moment expansion is performed by derivation with respect to the dust spectral index β. Here, we present the results from an additional simulated data set in order to shed light on the capacity of the moments expansion on β to capture the complexity arising specifically from the spatial variations of the dust temperature $T(\widehat{\mathit{n}})$. Our SIM1 data set has constant $\beta (\widehat{\mathit{n}})={\beta }_0$ and $T(\widehat{\mathit{n}})=T_0$ over the sky, while our SIM3 data set has spatial variations of both these parameters; here we introduce SIM6, a data set simulated with a varying spectral index $\beta (\widehat{\mathit{n}})$ and a constant temperature T₀:

SIM6: $M_{{\nu }_i}^{\text{SIM6},\eta }=N_{{\nu }_i}^{\eta }+S_{{\nu }_i}^{{\text{D}}_4}$,

where S^D₄ is the dust model with varying spectral index $\beta (\widehat{\mathit{n}})$ and a constant temperature T₀, presented in Sect. 4.1.2. The SIM6 simulated data set underwent the same fits as for the other sets presented in the main body of this article.

Figure E.1 shows the reduced χ² results of the MBB and the dust moment fits up to third order for the SIM6 simulations. They globally show similar trends as the ones computed for SIM3 (see Fig. 1). The reduced χ² values for the MBB are significantly higher than that of SIM3, while higher order fits give slightly better χ².

Fig. E.1. Reduced χ2 of the fit of Eq. (5) at zero (solid black), first (dashed yellow), second (dashed-dotted blue), and third order (dotted red) for the SIM6 ($\beta (\widehat{\mathit{n}}),T_0$, orange crosses).

Pushing further the comparison between SIM3 and SIM6, we can see from Fig. E.2 that most of the parameters from the fits are similar, except for ${\beta }_0^{\text{corr}}$ and ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$. SIM6 shows significantly higher values than SIM3 for ${\beta }_0^{\text{corr}}$ and ${\mathcal{M}}_{\mathcal{l}}^{A_{\text{D}}{\omega }_2}$.

Fig. E.2. Comparison between SIM3 (red diamonds) and SIM6 (orange crosses) fitted parameters for the corrected spectral index ${\beta }_0^{\text{corr}}$ and the higher order moments (${\mathcal{M}}_{\mathcal{l}}^{ab}$ (143 × 545)∈{AD, ω1, ω2, ω3}) as a function of the multipole ℓ.