Issue 
A&A
Volume 647, March 2021



Article Number  A52  
Number of page(s)  16  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201937367  
Published online  08 March 2021 
Dust moments: towards a new modeling of the galactic dust emission for CMB Bmodes analysis
^{1}
IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France
email: jonathan.aumont@irap.omp.eu
^{2}
Laboratoire de Physique de l’ENS, ENS, Université PSL, CNRS, Sorbonne Université, Université ParisDiderot, Paris, France
^{3}
Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK
^{4}
School of Physical Sciences, National Institute of Science Education and Research, HBNI, Jatni 752050, Odisha, India
Received:
19
December
2019
Accepted:
4
December
2020
The characterization of the spectral energy distribution (SED) of dust emission has become a critical issue in the quest for primordial Bmodes. The dust SED is often approximated by a modified black body (MBB) emission law but the extent to which this is accurate is unclear. This paper addresses this question, expanding the dust SED at the power spectrum level. The expansion is performed by means of moments around the MBB law, related to derivatives with respect to the dust spectral index. We present the mathematical formalism and apply it to simulations and Planck total intensity data, from 143 to 857 GHz, because no polarized data are yet available that provide the required sensitivity to perform this analysis. With simulations, we demonstrate the ability of highorder moments to account for spatial variations in MBB parameters. Neglecting these moments leads to poor fits and a bias in the recovered dust spectral index. We identify the main moments that are required to fit the Planck data. The comparison with simulations helps us to disentangle the respective contributions from dust and the cosmic infrared background to the highorder moments, but the simulations give an insufficient description of the actual Planck data. Extending our model to cosmic microwave background Bmode analyses within a simplified framework, we find that ignoring the dust SED distortions, or trying to model them with a single decorrelation parameter, could lead to biases that are larger than the targeted sensitivity for the next generation of CMB Bmode experiments.
Key words: cosmic background radiation / submillimeter: ISM / methods: data analysis
© A. Mangilli et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The precise characterization of the properties of the polarized dust emission from our Galaxy is a crucial issue in the quest for the primordial Bmodes 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 Bmode 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 Bmode 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 Balloonborne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) observations have shown that the polarization fraction is remarkably constant from the farinfrared 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 soughtafter CMB Bmode 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 crosscorrelation 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 multifrequency 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 adhoc, 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 Bmodes.
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 crosspower spectra, which are highly relevant to CMB Bmode 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 (PlanckHFI) 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, PlanckHFI 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 Bmodes 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 crosspower 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 Bmodes 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 crossspectra 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 crosspower 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 momentexpansion dust SED parameterizations in Sect. 2.1. In Sects. 2.2 and 2.3 we present the spherical harmonics and the crosspower 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 RayleighJeans 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 singletemperature 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 β_{0} and a constant temperature T_{0} across the sky. At a given frequency ν the dust intensity map takes the form:
where is in this case a frequencyindependent dust intensity template map, ν_{0} a reference frequency, and B_{ν}(T_{0}) is the Planck law for the temperature T_{0}. 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 RayleighJeans 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):
where now 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 at a frequency ν now reads:
where 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^{1}), and .
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:
where^{2}β_{0}(ℓ) refers to the effective value of the dust spectral index β_{0} 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 . We note that the spherical harmonics moments (ω_{i})_{ℓm} are not the same as the spatial moments , 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 lineofsight or beam averaging. We also stress that the lineofsight and beamaveraging effects, for realworld experiments, can never be avoided, such that should already be interpreted as an averaged dustamplitudeweighted quantity.
2.3. Dust SED parametrization in the crosspower spectra
Based on Eq. (4), we can compute the crossspectrum between the maps observed in the frequency bands ν_{1} and ν_{2}. Up to third order in terms of β derivatives, this takes the form:
where the moment crosspower spectra, , {a, b}∈{A_{D}, ω_{1}, ω_{2}, ω_{3}} are defined as^{3}:
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 crossspectrum, , as the moments normalized to the dust amplitude spectrum and rescaled by the corresponding frequencydependent numerical factors for the ν_{i} × ν_{j} crossspectrum, as:
where c^{ab}(ν_{i}, ν_{j}, ν_{0}) are the numerical coefficients which involve the sum and/or the product of the ln(ν_{i}/ν_{0}) and ln(ν_{j}/ν_{0}) terms, as appearing in Eq. (5). In this way, the functions of Eq. (7) show, for a given crossspectrum, the effective contribution of each moment to the departure of the standard MBB SED. We stress that, in Eq. (5), the dust spectral index β_{0}(ℓ) 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, ), Eq. (5) becomes the crosspower 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 crosspower 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 adhoc choice.
We also note that the normalized firstorder term can be interpreted as a correction to the MBB spectral index β_{0}(ℓ) needed to recover the “true” β(ℓ) when spatial and/or lineofsight variations are present. We therefore define:
where the Δβ_{0}(ℓ) 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 scaledependent 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 β_{0}(ℓ) in Eq. (8) is related to a powerspectrum weighted average of (which, strictlyspeaking, itself is a lineofsightaveraged quantity). A similar averaging process was discussed for the temperature power spectrum stemming from the relativistic SunyaevZeldovich effect (Remazeilles et al. 2019). The higher order moment terms in Eq. (5) quantify (3D) crosscorrelations 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 lineofsight variations of the dust spectral index β at the powerspectrum 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 crossspectra data vector 𝒟_{ℓ}, gathering all the cross spectra that can be computed from a set of maps M_{νi} at frequencies ν_{i}, i ∈ {1, …, N}:
The analysis is divided in three steps:
– Step 1: We fit for each multipole ℓ the two parameters β_{0}(ℓ) and (T_{0} is fixed). Zeroorder or MBB fits, in the following, refer to this first step
– Step 2: We fix β_{0}(ℓ) 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 , fix it and redo the fits as in step 2. In practice, step 3 may need to be repeated in order to ensure that Δβ_{0}(ℓ) and thus 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: , , , , , , , , , . This “threestep” 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 β_{0}(ℓ) due to SED averaging effects, and therefore providing a description of the scale dependence of the frequency decorrelation.
In practice, we perform a LevenbergMarquardt leastsquares minimization of r = y − m_{i}, where, at each multipole bin, y ≡ 𝒟_{ℓ} is the input crossspectra data vector of Eq. (9) and is the model vector:
For step 1, the model refers to the dust crossspectra function truncated at zero order, m_{i} ≡ m_{0}, that is, the MBB emission with constant temperature and spectral index β_{0}(ℓ). For steps 2 and 3, the model refers to the dust moment crossspectra function of Eq. (5), truncated at orders 1, 2, and 3.
In order to take into account the correlations between the crossspectra 𝒟^{νi × νj} and 𝒟^{νk × νl}, the crossspectra covariance matrix ℂ is included in the fit:
The computation of this crossspectra covariance matrix from simulations is described in detail, for our applications, in Appendix B. The reduced χ^{2} to be minimized is then defined as:
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 highestfrequency PlanckHFI 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 freefree 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^{4} pixelization (Górski et al. 2005).
4.1. Simulated map data sets
For each simulation type, the frequency channel maps 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 EndtoEnd noise maps (, 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^{5}).
4.1.2. Dust component
The dust component is built from a dust intensity map template at 353 GHz, defined, in MJy sr^{−1} units, as:
where B_{ν = 353 GHz}(T_{0} = 19.6 K) is the Planck function at 353 GHz for the temperature T_{0} = 19.6 K in W m^{2} sr^{−1}, and is the dust optical depth map at 353 GHz. For all the simulations, we use the 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)^{6}. The GNILC analysis produced allsky 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 at 353 GHz as defined in Eq. (13) is shown in the top panel of Fig. B.2.
The coefficients used to rescale the dust template at 353 GHz to a different frequency ν_{i} are defined as:
meaning that the dust map at the frequency ν_{i} is:
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. : simulations with constant dust temperature ( K) and spectral index ().
2. : simulations with Gaussian variations of the dust spectral index and fixed temperature ( K). The Gaussian spectral index map is a single random realization (i.e., the same realization for all the simulations) of the normal distribution with a mean β_{0} = 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 map is shown in the top panel of Fig. B.3.
3. : simulations using Planck sky maps of the dust spectral index and temperature. The dust spectral index map and the temperature map are those derived from the MBB fit of the Planck GNILC dust total intensity maps (Planck Collaboration Int. XLVIII 2016). For illustration, the map is shown in the bottom panel of Fig. B.3.
4. : simulations using the GNILC Planck sky map of the dust spectral index and a constant dust temperature T_{0} = 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 PlanckHFI reference frequencies. For comparison, the PR3 data will be colorcorrected to account for the PlanckHFI 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 . For this, we use the multifrequency CIB simulation of the Planck Sky Model (PSM) version 1.9 described in Delabrouille et al. (2013)^{7}. 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: ,

(b)
SIM2: ,

(c)
SIM3: ,

(d)
SIM4: .
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 PlanckHFI frequency channel maps we use in this work. As the SMICA map resolution is different from that of the PlanckHFI frequency channel maps, we first deconvolve the SMICA CMB map from its beam function and then convolve it with the beam of individual PlanckHFI frequency channel maps before subtraction. All the required information to do so is available in the PLA^{5}.
4.3. Crosspower spectra computation
The crossangular 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 crossspectra we use the Xpol code described in Tristram et al. (2005). The Xpol code is a pseudoC_{ℓ} 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 crossspectra^{8} between the five PlanckHFI channels from 143 to 857 GHz, as described in Sect. 3, and store them in the data vector 𝒟_{ℓ} (see Eq. (9)). The crossspectra vector is binned in 15 bins of multipoles^{9} of size Δℓ = 20 in the range most relevant for CMB primordial Bmodes analysis (ℓ ∈ {20, 300}), such that 𝒟_{b} ≡ ⟨𝒟_{ℓ}⟩_{ℓ∈b}. As we work only with the binned version of the crossspectra, 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 autocorrelation bias and to reduce the level of correlated systematic errors, we compute the crossspectra from data split maps (the Planck halfmission maps, HM). We explain how we combine halfmission maps and compute the covariance matrix of the crossspectra in Appendix B.
After the computation of the crossspectra 𝒟_{ℓ}, we subtract, from every data set (PR3 and simulations), the averaged crossspectrum computed from the 300 Planck EndtoEnd 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 colorcorrection to the PR3 data set crossspectra to get rid of PlanckHFIspecific calibration effects as described in Planck Collaboration IX (2014). The units of the total intensity crossspectra for all the data sets are [(MJy sr^{−1})^{2}].
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 threestep 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 in Sect. 5.1, on the dust spectral index β_{0}(ℓ) and its leading order correction Δβ_{0}(ℓ) 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 χ^{2}(ℓ) plots, one per data set, shown in Fig. 1. In each plot, we show the reduced χ^{2}(ℓ) 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 (dasheddotted 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 χ^{2} 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 thirdorder fit is needed for the χ^{2}(ℓ) to reach unity. For SIM2, fitting with higher order terms is required to get a reduced χ^{2}(ℓ)∼1. First order gives a fair χ^{2}(ℓ) (except at very low ℓ) and second order is needed to get a χ^{2}(ℓ) 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 signaltonoise ratio. For SIM3, a moment expansion up to secondorder is required to get a fair fit and up to third order to reach χ^{2}(ℓ)∼1, while for the SIM4 set that includes the CIB, the third order is needed to have a fair reduced χ^{2}(ℓ). For the PR3 data, the reduced χ^{2}(ℓ) 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 for the thirdorder 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 dustonly 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 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 β_{0}(ℓ) derived from the MBB fit, and the bottom plot the corrected dust spectral index 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: from step 3 of the model fit. The symbols for the different data sets are as in Fig. 2. 
We find that β_{0}(ℓ) matches the input spectral index of the simulation β_{0} = 1.59 for the SIM1 set. For SIM2, β_{0}(ℓ) values are slightly larger than the input value at ℓ ≳ 100 but the small difference is corrected when computing .
Our results are less straightforward for the additional sets. For SIM3, we find with no systematic dependence with ℓ. For comparison, we computed the spectral index for a MBB from the ratio between the crossspectra 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_{0} = 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 .
The comparison of SIM3 and SIM4 in Fig. 3 shows that the CIB lowers both β_{0}(ℓ) and for ℓ > 100. Above this multipole, the difference between SIM3 and SIM4 values of increases steadily with ℓ, as does the CIB contribution. The values of 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 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 move in the opposite direction for SIM4 and PR3. While the two sets of values roughly match for the secondorder 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 crossspectra 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 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. Highorder moments
In this section, we discuss the results for the highorder (first to thirdorder) moments in Eq. (5). We choose to show the moment functions defined in Eq. (7) with ν_{i} = 143 GHz and ν_{j} = 545 GHz. We stress that this specific choice of frequencies only impacts the scaling prefactor c^{ab}(ν_{i}, ν_{j}, ν_{0}) 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 firstorder moment function , which is made compatible with zero through iteration on Δβ_{0}(ℓ) (see Sect. 3).
Fig. 4.
Step 3 moment functions as a function of the multipole ℓ defined in Eq. (7) for the 143 × 545 crossspectrum. The plots refer, from top left to bottom right, to , , , , , , , , and . 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 highorder moments. For the SIM2 data set, the function is detected with an amplitude increasing with ℓ, while the other moments are consistent with zero. For SIM3, we find that the moment is nearly scaleindependent and is significantly larger than that of SIM2 at low ℓ. For this set, two additional moment functions are detected, namely and more marginally . 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 . 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 thirdorder) are detected with an absolute amplitude larger than that measured on the simulated data sets. The moment function shows a similar amplitude as for SIM3 and SIM4 for ℓ ≳ 100, while at larger angular scales it does not. For the and , the PR3 data set has an overall ℓ behavior close to that of SIM4 but with an increased absolute value. The , , and of PR3 match those of SIM4 for ℓ ≲ 70 but progressively deviate from them for higher multipoles. Finally, we point out that the 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, highorder 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 highorder moments, we obtain a significantly different value that cancels the firstorder moment Aω_{1}.
– 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 , and . 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 highorder 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 highorder 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 tensortoscalar ratio measurement
We finally discuss the potential impact of our results on the measurement of the tensortoscalar 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 Bmodes 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 χ^{2} 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 Bmodes.
We consider the SIM3 crossspectra fit with Eq. (5) at different orders in the moment expansion. We look here at the relative difference between the SIM3 data set crossspectra and the fitted model of Eq. (5), for every crossspectra between frequencies ν_{i} and ν_{j}. We focus on the multipole bin centered at ℓ_{0} = 80, as this scale corresponds to the CMB primordial Bmodes peak. This relative difference reads:
This relative difference is displayed in Fig. 5. Focusing on the 143 × 143 crossspectrum, which is a frequency channel indicative of typical CMB Bmode experiments, the simple MBB fit leaves a Δ𝒟_{ℓ0}(143 × 143) = 10.9% residual, the firstorder fit leaves Δ𝒟_{ℓ0}(143 × 143) = 3.2%, the secondorder fit leaves Δ𝒟_{ℓ0}(143 × 143) = 0.5% and the thirdorder fit leaves Δ𝒟_{ℓ0}(143 × 143) = 0.06%.
Fig. 5.
Upper panel: mean SED in MJy^{2} sr^{−1} for the 100 SIM3 crossspectra as a function of the effective frequency (, in GHz) for the multipole bin centered at ℓ = 80. The MBB (solid black), first (dashed yellow), second (dasheddotted blue), and thirdorder (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 thirdorder (red diamonds) best fits. 
Let us now consider a future CMB Bmode experiment that has the same frequency channels and a S/N for the dust Bmodes that is similar to those of PlanckHFI for the dust intensity (e.g., LiteBIRD, Hazumi et al. 2019). We know from Planck Collaboration Int. XXX (2016) that the dust Bmode power spectrum computed on the LR42 mask has an amplitude of A_{D}(ℓ_{0}) = 78.6 μK^{2} at 353 GHz. If we convert this amplitude into the requivalent 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 Bmode experiment looking at half the sky could find an analysis bias of Δr = 0.109 ⋅ r_{D} = 0.20 by assuming that the dust Bmodes follow a MBB SED, Δr = 0.032 ⋅ r_{D} = 0.06 by assuming an firstorder moment expansion, Δr = 0.005 ⋅ r_{D} = 0.009 by assuming a secondorder moment expansion, and Δr = 0.0006 ⋅ r_{D} = 0.001 by assuming a thirdorder 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 Bmodes would lead to Δr = 0.109 ⋅ r_{D} = 0.01, a decorrelation or firstorder analysis would lead to Δr = 0.032 ⋅ r_{D} = 4 × 10^{−3}, a secondorder analysis to Δr = 0.005 ⋅ r_{D} = 5 × 10^{−4}, and a thirdorder analysis to Δr = 0.0006 ⋅ r_{D} = 7 × 10^{−5}.
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 Bmode analyses, targeting r values down to 10^{−3} and beyond.
Estimates of the tensortoscalar ratio bias Δr at 150 GHz for the LR42 and the BICEP2/Keck regions, in the case of the SIM3 crossspectra SED fitted assuming a MBB and a first, second, and or thirdorder SED moment expansion.
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 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^{−3}, we need to consider all moments with an absolute amplitude larger than 10^{−3}. 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, dustdominated angular scales of the intensity PR3 data set additionally stress the need for a thirdorder expansion of the polarized dust SED in order to reach the accuracy targeted by future CMB Bmode experiments.
7. Summary and conclusions
In this paper, we present a model that describes the Galactic dust SED for total intensity at PlanckHFI 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 highorder 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 crossspectra computed from the combination of CMBcorrected 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 Bmode 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 adhoc parameter, could lead to biases larger than the accuracy of the component separation required to search primordial Bmodes down to a tensortoscalar ratio r = 10^{−3}.
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 Bmode 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 dustdedicated frequency channels are needed in order to perform second and thirdorder moment expansion fits, respectively.
Additional difficulties for Bmode searches could arise from changes in polarization angles across frequencies, which would make the decomposition of polarized dust emission in E and Bmodes frequencydependent. 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 Bmodes. 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: ANR17CE310022). A.M. acknowledges the support of the DIAASTRO 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
 Abitbol, M. H., Chluba, J., Hill, J. C., & Johnson, B. R. 2017, MNRAS, 471, 1126 [Google Scholar]
 Ashton, P. C., Ade, P. A. R., Angilè, F. E., et al. 2018, ApJ, 857, 10 [Google Scholar]
 Chiang, Y.K., & Ménard, B. 2019, ApJ, 870, 120 [Google Scholar]
 Chluba, J. 2016, MNRAS, 460, 227 [Google Scholar]
 Chluba, J., & Sunyaev, R. A. 2012, MNRAS, 419, 1294 [NASA ADS] [CrossRef] [Google Scholar]
 Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195 [Google Scholar]
 Chluba, J., Abitbol, M. H., Aghanim, N., et al. 2019, ArXiv eprints [arXiv:1909.01593] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J. B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desjacques, V., Chluba, J., Silk, J., de Bernardis, F., & Doré, O. 2015, MNRAS, 451, 4460 [Google Scholar]
 Draine, B. T., & Hensley, B. 2013, ApJ, 765, 159 [Google Scholar]
 Gandilo, N. N., Ade, P. A. R., Angilè, F. E., et al. 2016, ApJ, 824, 84 [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hanany, S., Alvarez, M., Artis, E., et al. 2019, ArXiv eprints [arXiv:1902.10541] [Google Scholar]
 Hazumi, M., Ade, P. A. R., Akiba, Y., et al. 2019, J. Low Temp. Phys., 194, 443 [Google Scholar]
 Hensley, B. S., & Bull, P. 2018, ApJ, 853, 127 [Google Scholar]
 Juvela, M., Montillaud, J., Ysard, N., & Lunttila, T. 2013, A&A, 556, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. Lett., 78, 2058 [Google Scholar]
 Keck Array & BICEP2 Collaborations (Ade, P. A. R. et al.) 2018, Phys. Rev. Lett., 121, 221301 [CrossRef] [PubMed] [Google Scholar]
 Kogut, A., Abitbol, M. H., Chluba, J., et al. 2019, BAAS, 51, 113 [Google Scholar]
 Krachmalnicoff, N., Baccigalupi, C., Aumont, J., Bersanelli, M., & Mennella, A. 2016, A&A, 588, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XVII. 2014, A&A, 566, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXII. 2015, A&A, 576, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. L. 2017, A&A, 599, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poh, J., & Dodelson, S. 2017, Phys. Rev. D, 95, 103511 [Google Scholar]
 Polnarev, A. 1985, Sov. Astron., 29, 607 [NASA ADS] [Google Scholar]
 Remazeilles, M., Bolliet, B., Rotti, A., & Chluba, J. 2019, MNRAS, 483, 3459 [Google Scholar]
 Rotti, A., & Chluba, J. 2021, MNRAS, 500, 976 [Google Scholar]
 Sathyanarayana Rao, M., Subrahmanyan, R., Udaya Shankar, N., & Chluba, J. 2015, ApJ, 810, 3 [Google Scholar]
 Seljak, U., & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054 [Google Scholar]
 Shariff, J. A., Ade, P. A. R., Angilè, F. E., et al. 2019, ApJ, 872, 197 [Google Scholar]
 Sheehy, C., & Slosar, A. 2018, Phys. Rev. D, 97, 043522 [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1970a, Comm. Astrophys. Space Phys., 2, 66 [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1970b, Ap&SS, 7, 20 [Google Scholar]
 Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90 [Google Scholar]
 Tristram, M., MaciasPerez, J., Renault, C., & Santos, D. 2005, MNRAS, 358, 833 [Google Scholar]
 Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zeldovich, Y. B., & Sunyaev, R. A. 1969, Ap&SS, 4, 301 [NASA ADS] [CrossRef] [Google Scholar]
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 multicomponent simulations. We provide an illustrative example of the foreground templates at 353 GHz in MJy sr^{−1} 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 maps used in the simulations. The top panel of Fig. B.3 shows the map used to generate the SIM2 dust simulations with Gaussian variations with Δβ = 0.1 around β_{0} = 1.59. The bottom panel of Fig. B.3 shows the GNILC map used to generate SIM3 simulations (Planck Collaboration Int. XLVIII 2016).
Appendix B: Crossspectra 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 crosspower 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 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.
Crossspectra 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 crossspectra as defined in Eq. (B.1). Middle panel: correlation matrix computed for our toymodel of the correlation presented in Appendix B.2. Bottom panel: correlation matrix computed from the crossspectra as defined in Eq. (B.6). 
B.1. Definitions
In order to avoid the noise autocorrelation bias and to reduce the level of correlated systematic errors, we compute the crossspectra from data split maps (the 2 Planck halfmission maps, HM, namely HM1 and HM2). The general philosophy would be to reproduce Planck Collaboration Int. XXX (2016) to construct the crosspower spectra from two frequency maps M_{νi} and M_{νj}:
The idea of doing the sum of the four crossspectra 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 crossspectra, as we see in the following. To preserve the statistical independence between the crossspectra, 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 halfmission 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 crossspectra of our data sets, we build the correlation matrix from the covariance matrix of Eq. (11):
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 crossspectra. This highlights that the crossspectra, 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 crossspectra in Eq. (B.1) to minimize these correlations.
B.2. Expected correlations from a toy model
In a highS/N mixture of interstellar dust and instrumental noise, the correlation between frequency channels crossspectra 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:
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 highS/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^{10}:
Let us now consider as an example the specific correlation between 2 Planck crossspectra, 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_{217} × N_{143}, D_{353} × N_{143}) 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:
For a typical dust MBB SED with β_{0} = 1.59 and T_{0} = 19.6 K and the PlanckHFI 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 crossspectra. 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 crossspectra covariances in the correlation matrix ℝ, we can build the full toymodel 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 toymodel 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 crossspectra of our data sets, we change the definition of Eq. (B.1) in the case i ≠ j:
This is the definition of the data sets crossspectra we adopt in this paper. The correlation matrix built for the SIM1 data set from this definition of the crossspectra 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 crossspectra 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 crossspectra is N_{d.o.f.} = 5. When defining the crossspectra as in Eq. (B.6), it becomes N_{d.o.f.} = 11, allowing us to fit the SED moment expansion of the crossspectra 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 thirdorder dust SED moment expansion. Here, we present the corrected dust amplitude , the corrected spectral index , and the relevant moment functions ℳ_{ℓ} in the case where Eq. (5) is truncated at first and second order in our threestep fitting procedure.
These other truncating order results are displayed in Figs. C.1–C.6 for , , , , and , respectively.
Fig. C.1.
Corrected dust amplitude 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). 
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, 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: ,
where 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 subdominant with respect to the CIB by at least three orders of magnitude.
Figure D.1 shows the reduced χ^{2} 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 , the dust spectral index β_{0}(ℓ), the corrected spectral index , and the dust moments functions 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 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 (dasheddotted 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} = (X^{SIM5} − X^{SIM4})/X^{SIM4} with 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 PlanckHFI 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 . Our SIM1 data set has constant and 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 and a constant temperature T_{0}:
SIM6: ,
where S^{D4} is the dust model with varying spectral index and a constant temperature T_{0}, 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 χ^{2} 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 χ^{2} values for the MBB are significantly higher than that of SIM3, while higher order fits give slightly better χ^{2}.
Fig. E.1.
Reduced χ^{2} of the fit of Eq. (5) at zero (solid black), first (dashed yellow), second (dasheddotted blue), and third order (dotted red) for the SIM6 (, 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 and . SIM6 shows significantly higher values than SIM3 for and .
Fig. E.2.
Comparison between SIM3 (red diamonds) and SIM6 (orange crosses) fitted parameters for the corrected spectral index and the higher order moments ( (143 × 545)∈{A_{D}, ω_{1}, ω_{2}, ω_{3}}) as a function of the multipole ℓ. 
The known dust emission β–T anticorrelation (e.g., Juvela et al. 2013) could be an explanation for the differences between results on SIM6 and SIM3. SIM6, with a constant T_{0} on the sky, shows a larger variability in SED due to the β spatial variations, whereas in SIM6 the spatial variations of T – anticorrelated to those of β – tend to compensate the SED variability, hence the larger χ^{2} when fitting the MBB.
All Tables
Estimates of the tensortoscalar ratio bias Δr at 150 GHz for the LR42 and the BICEP2/Keck regions, in the case of the SIM3 crossspectra SED fitted assuming a MBB and a first, second, and or thirdorder SED moment expansion.
All Figures
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 (dasheddotted blue), and third order (dotted red) are displayed. 

In the text 
Fig. 2.
Amplitude of the dust spectrum 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. 

In the text 
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: from step 3 of the model fit. The symbols for the different data sets are as in Fig. 2. 

In the text 
Fig. 4.
Step 3 moment functions as a function of the multipole ℓ defined in Eq. (7) for the 143 × 545 crossspectrum. The plots refer, from top left to bottom right, to , , , , , , , , and . The symbols refer to the different data sets: SIM1 (green circles), SIM2 (yellow squares), SIM3 (red diamond), SIM4 (blue stars), and PR3 (black triangles). 

In the text 
Fig. 5.
Upper panel: mean SED in MJy^{2} sr^{−1} for the 100 SIM3 crossspectra as a function of the effective frequency (, in GHz) for the multipole bin centered at ℓ = 80. The MBB (solid black), first (dashed yellow), second (dasheddotted blue), and thirdorder (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 thirdorder (red diamonds) best fits. 

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

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

In the text 
Fig. B.3.
Dust spectral index map 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). 

In the text 
Fig. B.4.
Crossspectra 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 crossspectra as defined in Eq. (B.1). Middle panel: correlation matrix computed for our toymodel of the correlation presented in Appendix B.2. Bottom panel: correlation matrix computed from the crossspectra as defined in Eq. (B.6). 

In the text 
Fig. C.1.
Corrected dust amplitude 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). 

In the text 
Fig. C.2.
Same as Fig. C.1 but for . 

In the text 
Fig. C.3.
Same as Fig. C.1 but for , normalized for the 143 × 545 GHz crossspectrum. 

In the text 
Fig. C.4.
Same as Fig. C.1 but for , normalized for the 143 × 545 GHz crossspectrum. 

In the text 
Fig. C.5.
Same as Fig. C.1 but for , normalized for the 143 × 545 GHz crossspectrum. 

In the text 
Fig. C.6.
Same as Fig. C.1 but for , normalized for the 143 × 545 GHz crossspectrum. 

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

In the text 
Fig. D.2.
Relative difference between SIM5 and SIM4 fitted parameters. Δ_{X} = (X^{SIM5} − X^{SIM4})/X^{SIM4} with as a function of the multipole ℓ. 

In the text 
Fig. E.1.
Reduced χ^{2} of the fit of Eq. (5) at zero (solid black), first (dashed yellow), second (dasheddotted blue), and third order (dotted red) for the SIM6 (, orange crosses). 

In the text 
Fig. E.2.
Comparison between SIM3 (red diamonds) and SIM6 (orange crosses) fitted parameters for the corrected spectral index and the higher order moments ( (143 × 545)∈{A_{D}, ω_{1}, ω_{2}, ω_{3}}) as a function of the multipole ℓ. 

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