Dust moments: towards a new modelling of the galactic dust emission for CMB B-modes analysis

The characterization of the spectral energy distribution (SED) of dust emission has become a critical issue in the quest for primordial B-modes. The dust SED is often approximated by a modified blackbody (MBB) emission law but to what extent is it accurate? This paper addresses this question expanding the dust SED, at the power spectrum level, in 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 today no polarized data provide the required sensitivity to perform this analysis. With simulations, we demonstrate the ability of high order moments to account for spatial variations in MBB parameters. Neglecting these moment 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 help us to disentangle the respective contributions from dust and the Cosmic Infrared Background to the high order moments, but the simulations give an insufficient description of the actual Planck data. Extending our model to CMB B-mode 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 larger than the targeted sensitivity for the next generation of CMB B-mode experiments.


Introduction
The precise characterization of the properties of the polarized dust emission from our Galaxy is a crucial issue for the quest of the Cosmic Microwave Background (CMB) primordial Bmodes. If measured, this faint cosmological signal imprinted by the primordial gravitational wave background, would be an evidence of the inflation epoch and 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, an accurate determination of diffuse CMB B-mode foregrounds, among which the polarized Galactic dust emission -dominating at observing frequencies 70 GHz (see e.g. Krachmalnicoff et al. 2016;Planck Collaboration et al. 2018a) -is required to get an unbiased estimate of the ratio r between tensor and scalar primordial perturbations, a parameter of an 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 has shown that the SED of the dust emission for total intensity and polarization can be fitted by a modified blackbody law, hereafter MBB (Planck Collaboration et al. 2014b,c, 2015, 2018a. Maps of the dust MBB spectral indices and temperatures have been fitted to the total intensity Planck data (see e.g. Planck Collaboration et al. 2014aCollaboration et al. , 2016a. 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 (SNR).
Send offprint requests to: jonathan.aumont@irap.omp.eu Analyzing total intensity data is however 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 millimetre wavelengths, suggesting that the polarized and total dust emission arise predominantly from a single grain population (Planck Collaboration et al. 2015;Gandilo et al. 2016;Ashton et al. 2017;Guillet et al. 2018;Shariff et al. 2019;Planck Collaboration et al. 2018a).
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 overcome. 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 et al. 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 maps 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 et al. 2017) has been later dismissed (Sheehy & Slosar 2018). This analysis was further extended by Planck Collaboration et al. (2018a), performing a global multi-frequency fit of the polarized Planck HFI channels with a spectral model that in-cludes decorrelation. They 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 (BICEP2 Collaboration et al. 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 dust polarized emission has thus become a main challenge 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 spherical harmonics expansion. 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 non-linear 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 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, getting rid of noise bias and uncorrelated systematic effects. Similar extensions could also play an important role for 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 SNR 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. Sect. 3 details our methodology and implementation of the dust moments analysis. We present the simulations and the actual Planck data on which we fit our spectral model in Sect. 4; the results of the fits are 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 four 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.

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 powerful to account for SED distortions arising from the various averaging effects. We first recall the usual and the moment-expansion dust SED parametrizations in Sect. 2.1. In Sect. 2.2 and Sect. 2.3 we present the spherical harmonics and the cross-power spectra moment expansion that will be 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 et al. 2016a). This can be easily generalized to include derivatives with respect to the dust temperature . However, in the following analysis, we only use the spectral index moment expansion, as temperature and βvariations can have similar effects on the dust SED built from microwave and sub-millimetre data, in the Rayleigh-Jeans regime.

Modified blackbody (MBB) formalism
The commonly used parametrization for a single-temperature dust spectrum is a modified blackbody (MBB) emission law. We first consider the MBB parametrization without any spatial variations of the SED, i.e. with a constant spectral index β 0 and a constant temperatude T 0 across the sky. At a given frequency ν the dust intensity map I D (ν,n) takes the form: where A D (n) is in this case a frequency independent 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 3 dimensions of space. The first effect is that the effective spectral index β or temperature T in Eq.
(1) can vary spatially. 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 was accurately described by Eq. (1), the MBB no longer accurately describes the effective emission law.

MBB 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): clarity, in writing this equation, we ignore the variations along the line of sight.

MBB spectral index moment expansion
A more general and powerful parametrization of the dust SED has been 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 D (ν,n) at a frequency ν now reads: where ω i (n) = A D (n)∆β(n) i is the i th order moment map associated to the i th derivative of the MBB with respect to β (here in an expansion up to the 3 rd order 1 ), and ∆β(n) = β(n) − β 0 .

Dust SED parametrization in spherical harmonics
By definition, the decomposition of the dust map into spherical harmonics coefficients implies averages over various linesof-sight over the sky, which is mathematically equivalent to averaging pixels with different SEDs along the line of sight or in an instrumental beam . Therefore, we can use the spectral moment decomposition described in Eq. (3). It 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 will see in the following. The variables (ω i ) m (i ∈ {1, . . . , N}) are the spherical harmonic coefficient of the i th order moment map ω i (n). Note that the spherical harmonics moments (ω i ) m are not the same as the spatial moments ω i (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 already β(n) should be interpreted as an averaged dust-amplitude-weighted quantity.

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 ν 1 and ν 2 . Up to the third order in terms of β derivatives, it takes the form: D ν 1 ×ν 2 = I ν 1 (T 0 , β 0 ( ))I ν 2 (T 0 , β 0 ( )) I 2 ν 0 (T 0 , β 0 ( )) where the moment cross-power spectra, D ab , {a, b} ∈ {A D , ω 1 , ω 2 , ω 3 } are defined as 3 : In Eq. (5), we grouped terms according to the maximal derivative order in β that appears. This 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 3 terms, truncating after the second moment order means including the first 6 terms and so on. Since it will be useful in the following, we define here the moment functions for a given cross-spectrum, M ab (ν i , ν j ), as the moments D ab normalized to the dust amplitude spectrum D A D A D and re-scaled by the corresponding frequency-dependent numerical factors for the ν i × ν j cross-spectrum, 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 M 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 β 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 modelled. At the order 0 (keeping only the first element of the sum, D A D A 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, so that Eq. (5) provides us with a consistent and robust description of the dust SED distortions with respect to the modified blackbody 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 gives an efficient framework 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 et al. 2017;Vansyngel et al. 2017;BICEP2 Collaboration et al. 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 β 0 ( ) needed to recover the "true" β( ) when spatial and/or line-of-sight variations are present. We thus define: 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 scale dependent bias that arises from neglecting SED corrections. Adding the moment expansion allows us to eliminate the bias while having non-zero 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 power-spectrum weighted average of β(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 (3 dimensional) cross-correlations between the spectral indices . 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 physicallymotivated model for dealing with spatial and line-of-sight variations of the dust spectral index β at the power-spectrum level.

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 D , 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 3 steps: - Step 1: we fit for each multipole the two parameters β 0 ( ) and D A D A D (T 0 is fixed). Order 0 or MBB fits, in the following, refer to this first step -Step 2: we fix β 0 ( ) and then fit Eq. (5): 3 parameters are fitted at the 1 st order, 6 at the 2 nd order and 10 at the 3 rd order -Step 3: we update the value of the spectral index to β corr 0 ( ) = β 0 ( ) + ∆β 0 ( ), 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 D A D ω 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, will be: This "3-step" method for the dust moments fit allows to consistently search for spectral departures from the standard MBB quantifying, at each multipole, the corrections to the β 0 ( ) due to SED averaging effects, therefore providing a description of the scale dependence of the frequency decorrelation.
In practice, we perform a Levenberg-Marquardt leastsquares minimization of r = y − m i , where, at each multipole bin, y ≡ D is the input cross-spectra data vector of Eq. (9) and m i ≡ D D is the model vector: For the step 1, the model refers to the dust cross-spectra function truncated at the 0 th order, m i ≡ m 0 , that is, the modified blackbody emission with constant temperature and spectral index β 0 ( ). For step 2 and step 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 D ν i ×ν j and D ν k ×ν l , the cross-spectra covariance matrix C is included in the fit: The computation of this cross-spectra 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 dof is the number of degrees of freedom.

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, however 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 et al. 2014c).
In the following analysis, we will 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 (labelled SIM1, SIM2, SIM3 and SIM4) and the actual Planck intensity data (labelled PR3).

Simulated map data sets
For each simulation type, the frequency channel maps M SIM ν i 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, 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.

Noise component
The noise component N is the same for each simulation type. We use 300 Planck End-to-End noise maps (N η ν i , η ∈ {1, 300}, 300 realizations per frequency band ν i ) obtained from the FFP10 Planck simulations (which include the contributions of both the noise and the residual systematics). These maps are publicly available in the Planck Legacy Archive (PLA 4 ).

Dust component
The dust component is built from a dust intensity map template at 353 GHz, defined, in MJy sr −1 units, as: The coefficients used to re-scale the dust template at 353 GHz to a different frequency ν i are defined as: so that the dust map at the frequency ν i is: We define three types of dust simulation, with different level of complexity in the frequency scaling of the dust intensity in Eq. (14) and Eq. (15): 1. S D 1 ν i : simulations with constant dust temperature (T (n) = T 0 = 19.6 K) and spectral index (β(n) = β 0 = 1.59). 2. S D 2 ν i : simulations with Gaussian variations of the dust spectral index and fixed temperature (T (n) = T 0 = 19.6 K). The Gaussian spectral index map β(n) = N(β 0 , ∆β 2 ) 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 match roughly observed dispersion of the spectral index variations in the Planck data (Planck Collaboration et al. 2015, 2016c. The corresponding β(n) map is shown in the top panel of 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 et al. 2014d).

CIB and synchrotron
In one of our simulation sets, detailed below, we include a CIB component S CIB ν i . For this, we use the multi-frequency CIB simulation of the Planck Sky Model (PSM) version 1.9, described in Delabrouille et al. (2013) 6 . This CIB map is a random Gaussian realization matching the Planck measured CIB power spectra of Planck Collaboration et al. (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 here the simulations of this component, but comment on it in Appendix D.

Simulation types
From the components above, we produce 4 batches of 100 simulated intensity maps at 143, 217, 353, 545 and 857 GHz, numbered by the superscript η ∈ {1, 100}: 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.

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 will 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 et al. 2018b) 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 .

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 et al. (2016b) and used in several Planck publications. It leaves 42% of the sky for the analysis. The LR42 mask is apodized and includes a galactic mask, a point source mask and a CO mask, as described in Planck Collaboration et al. (2016b). It 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, correcting 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 7 between the five Planck-HFI channels from 143 to 857 GHz, as described in Sect. 3, and store them in the data vector D (see Eq. (9)). The cross-spectra vector is binned in 15 bins of multipoles 8 of size ∆ = 20 in the range most relevant for CMB primordial B-modes analysis ( ∈ {20, 300}), such that D b ≡ D ∈b . As we work only with binned version of the cross-spectra, in the rest of this document and for the sake of clarity, we do not make the distinction between b and , so that D refers to D b .
In order to avoid the noise auto-correlation bias and to reduce the level of correlated systematics, we compute the cross-spectra from data split maps (the Planck half-mission maps, HM). We explain how we combine half-mission maps and how we compute the covariance matrix of the cross-spectra in Appendix B.
After the computation of the cross-spectra D we subtract from every data set (PR3 and simulations), the averaged crossspectrum computed from the 300 Planck End-to-End simulations that include instrumental noise and systematic effects (see Sect. 4.1.1). This allows to correct for the small bias linked to residual systematics, 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 et al. (2014d). The units of the total intensity cross-spectra for all the data sets are [(MJy sr −1 ) 2 ].

Results for simulations and Planck data
In this section, we present the results of the fits for the five data sets (SIM1, SIM2, SIM3, SIM4 and PR3) described in Sect. 4. We fit the model D D on each data set D cross-spectra vector, independently in each multipole bin , using the 3-step method introduced in Sect. 3. up to order 1, 2 and 3 (for some of our data sets, step 3 need to be repeated once to ensure convergence). The maximum order of the expansion -3 -is set by the number of degrees of freedom in our data sets (see Appendix B).
The fit results are presented as follows. We describe the goodness of our fits and the dust amplitude spectrum D A D A D in Sect. 5.1, the dust spectral index β 0 ( ) and its leading order correction ∆β 0 ( ) in Sect. 5.2, and finally higher order moments fit in Sect. 5.3. In Sect. 5.4 we summarize and discuss these results.
For each set simulated data set, we present mean results averaged over 100 realizations, with error-bars corresponding to the standard deviation among realizations. For the PR3, we estimate error-bars propagating data uncertainties through the fits. These two approaches yield comparable error-bars.

Goodness of fit analysis and dust amplitude spectra
As the fits are performed independently in each multipole bin, the goodness of the fits is quantified with the reduced χ 2 ( ) plotted as a function of in Fig. 1. Each panel of the figure shows χ 2 ( ) for one data set fitted with the MBB law (step 1) and moments expansion in Eq. (5) truncated to the 1 st , 2 nd and 3 rd order (step 3). The reduced χ 2 quantifies the maximum order in Eq. (5), needed to describe each of the simulations and the PR3 data. This diagnostic allows us to compare simulations of increasing spectral complexity with the Planck PR3 data.
Obviously, the SIM1 set is well fitted by the MBB -order 0 of Eq. (5) -as in these simulations there are no variations of the MBB parameters. The fact that the reduced χ 2 of the MBB fit is close to 1 indicates that the correction for Planck residual systematics described in Sect. 4.3 is effective. There is a small hint for residual systematics in the lowest = 23 bin; indeed, for this bin a 3 rd order fit is needed for the χ 2 ( ) to reach unity. For the SIM2 set with random Gaussian variations of the SED, the MBB is not anymore a good fit and fitting with higher orders terms is required to get a reduced χ 2 ( ) of the order of 1. Order 1 gives a fair χ 2 ( ) (except at very low ) and order 2 is needed to get a χ 2 ( ) close to unity. For the SIM3 set a moment expansion up to order 2 is required to get a fair fit and order 3 to reach χ 2 ( ) ∼ 1. When including the CIB component in SIM4, the 3 rd order is needed to have a fair reduced χ 2 ( ). For the PR3 set, the reduced χ 2 ( ) is somewhat worst than those of the SIM4 set (except for the MBB fit which is slightly better), indicating that the Planck data have more spectral complexity than the SIM4 simulations.
The SIM3 and SIM4 simulations with dust SED variations based on Planck data and the PR3 data have huge χ 2 ( ) values for the MBB fits. We interpret these high values as due to the multipole averaging, which increases the signal to noise and, concomitantly, introduces spectral complexity. To our knowledge, the poorness of the dust intensity power spectrum fitting with a MBB was never reported before.
In Fig. 2, we present the dust amplitude spectra D A D A D from the step 3 fit at 3 rd order for the five data sets. For all data sets, these amplitudes are relatively insensitive to the order of the fit (see Fig. C.1) . This is expected, as they are not affected by the dust SED distortions, which are linked to variations of the dust spectral index and the dust temperature. The SIM1, SIM2 and SIM3 sets show indistinguishable results as they are built from the same dust spatial template and only differ in the modelling of the dust SED. The SIM4 and PR3 spectra are close to each other. Both depart from the dust-only simulations for multipoles 100, i.e. 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 in building the dust simulations are not fully free of CIB (Chiang & Ménard 2019).

Dust spectral index
The dust spectral indices derived from our fitting are presented for the five data sets in Fig. 3. The top panel displays the spectral indices β 0 ( ) derived from the MBB fit at the first step of our analysis, and the bottom panel the corrected dust spectral index β corr 0 ( ) ≡ β 0 ( ) + ∆β 0 ( ), which cancels the first order moment D Aω 1 when fitting the full model in Eq. (5) up to order 3. The correction depends on the truncation order of the fit as illustrated in Fig. C.2. For the SIM1 set, β 0 ( ) matches the input spectral index of the simulation β 0 = 1.59. The values for the SIM2 set are slightly larger than the input value at 100 but these small differences are corrected when computing β corr 0 ( ). The interpretation of our results is less straightforward for the SIM3, SIM4 and PR3 sets.
For the SIM3 set, the median value of the spectral index GNILC map (corrected to a mean temperature of 19.6 K) over the L42 mask is 1.6. The mean value derived from the ratio between the cross-spectra of the 217 and 353 GHz maps and the 353 GHz power spectrum, as in Planck Collaboration et al. (2018a), is lower: it ranges from 1.52 to 1.55 with no systematic dependence on . These latter values measured in harmonics space are close to those derived from our spectral analysis. The mean corrected dust spectral index β corr 0 ( ) SIM3 = 1.537 ± 0.003. The comparison of SIM3 and SIM4 in Fig. 3 shows that the CIB contribution lowers both β 0 ( ) and β corr 0 ( ) for > 100. Above this multipole, the difference between SIM3 and SIM4 values of β corr 0 ( ) increases steadily with . The -dependence of β corr 0 ( ) are remarkably similar for SIM4 and PR3 sets.  Fig. 3. Upper panel. Spectral index β 0 ( ) as a function of the multipole for the MBB fit (step 1). Bottom panel.
Step 3 corrected dust spectral index: β corr 0 ( ) = β 0 ( ) + ∆β 0 ( ). The symbols for the different data sets are as in Fig. 2. Its is interesting to compare our PR3 results with estimates presented in earlier determinations of the dust spectral index from Planck total intensity data. Planck Collaboration (2018) measured the dust spectral index, for a MBB SED, from the ratio between the 217 × 353 and the 353 × 353 cross-spectra, in 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 to 1.50, a bit larger than the mean value β corr 0 ( ) PR3 = 1.45 ± 0.01 we measure over the same multipole range. In Planck Collaboration (2018), the spectral index increases for higher multipoles up to 1.53 ± 0.01 in their = 140 − 169 bin. In our analysis, we observe an opposite trend with spectral index decreasing for > 100. These differences at < 100 reflect the impact of the spectral model on the determination of the dust spectral index. For our analysis, this is illustrated in Fig. C.2, which shows the dependence of β corr 0 ( ) with the truncation order of the model we fit.

High order moments
We finally present the results for the high (1 st to 3 rd ) order moments in Eq. (5). We choose to show the moment func-tions M ab (ν i , ν 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 , ν 0 ) in Eq. (7), and not the -dependence of the moments, nor the relative values between data sets for a given moment. The moment functions are presented, for the five data sets, each in one of the plots of Fig. 4. We present the moments from fits up to order 3 for all data sets, even when the reduced χ 2 drops to unity at a lower order (see Sect. 5.1). Nevertheless, we have checked that when this happens (e.g. for SIM1 and SIM2), the fit up to order 3 does not modify the lower order functions significantly (see Appendix C). The top left panel of Fig. 4 represents the first order moment function M A D ω 1 , that is made compatible with zero through iteration on ∆β 0 ( ) (see Sect. 3).
Even if the MBB law is a good fit for the SIM1 set, we performed the dust moments analysis up to the third order to demonstrate that the higher order moments are compatible with zero. As expected, we do not detect any moment function for SIM1, as shown in Fig. 4. This is a validation test of our method that allows us to conclude that residual systematics do not have a significant impact on high order moments. For the SIM2 data set, with Gaussian spatial variations of the dust SED, the M ω 1 ω 1 function is detected with an amplitude increasing with . The other moments are consistent with zero for this data set. The SIM3 data set, with realistic spatial variations of the dust SED, has a M ω 1 ω 1 moment close to scale independent and significantly larger than that of SIM2 at low . Furthermore, for this set two additional moment functions, the M ω 1 ω 2 and more marginally M A D ω 2 , are detected. Comparing moment functions for the SIM3 and SIM4 sets, we find that the CIB has a significant impact on several moments at 100, but for M ω 1 ω 1 . Almost all of the other moment functions increase with and are close to zero in the lowest -bins.
The analysis of the Planck PR3 data reveals more spectral complexity than observed for the SIM4 set. All the moment functions (from order 1 up to order 3) are detected with an absolute amplitude larger than that measured on the simulated data sets. The M ω 1 ω 1 moment function shows a similar amplitude as for SIM3 and SIM4 for 100, while at larger angular scales it does not. For the M A D ω 2 and M ω 2 ω 2 , the PR3 data set has an overall behavior close to that of SIM4 but with an increased absolute value. The M ω 1 ω 2 , M ω 1 ω 3 , M ω 2 ω 3 and M ω 3 ω 3 of PR3 are matching those of SIM4 for 70 but progressively deviate from them for higher multipoles. Finally, we point out that the M A D ω 3 moment function is the one with highest absolute amplitude for the PR3 set. Surprisingly it has the opposite sign, and a different -dependence, to the SIM4 moment.

Discussion
We list the main results of the fits before briefly discussing their interpretation.
-The goodness of the fit obtained for the simulations demonstrates the ability of the moments 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ω 1  -The comparison of the moments obtained for the SIM3 and SIM4 simulations, which only differ in the CIB component, provides insight on the CIB contribution. 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 order 3 -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 moments 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.
The moment functions in Fig. 4 are difficult to interpret at > 100 due to CIB contribution, but for lower multipoles, based on the comparison of spectra in Fig. 2, one could probably ignore the CIB and relate the results to dust emission properties. For the PR3 data at < 100, the three most significant moment functions in decreasing order are M A D ω 3 , M A D ω 2 and M ω 1 ω 1 . This result is unexpected and it is also interesting to point 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 compo-nents, which are not fully correlated on the sky (Draine & Hensley 2013;Guillet et al. 2018;Hensley & Bull 2018). However, it is not straightforward to provide a specific interpretation of these high order moments: their amplitude, scale dependence and hierarchy. One difficulty is that the moments expansion does not decompose the data in independent components: the high order moments depend on the expansion order as illustrated in Fig. C.3 to Fig. C.6 of Appendix C. Furthermore, the moments also quantify the SED averaging that occurs when going from pixel space to harmonic space. It is therefore hard 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 good match of the Planck data would be needed to quantify possible interpretations. However, putting together all higher order moment maps we can in principle reconstruct the probability distribution functions and maps of the dust spectral index, to be used for simulating the frequency decorrelation of dust maps.

Implications for the tensor-to-scalar ratio measurement
We finally give an insight on the potential impact of our results on the measurement of the tensor-to-scalar ratio r. We have 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 we observed and quantified to be absent in polarization. We here indeed assume that the dust B-modes SED has spectral departures from their mean MBB SED of the same relative order as the ones we observed in intensity and we derive the potential bias on r they would lead to by 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 level of departure from the MBB present in this simulated data set. One could have used SIM4 or PR3 data sets but, as we have seen in Sect. 5, it is not trivial to say 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 0 = 80, as this scale corresponds to the CMB primordial B-modes peak. This relative difference reads: This relative difference is displayed in Fig. 5. The upper panel presents the SIM3 cross-spectra SED at 0 = 80 for the 15 cross-spectra considered in the analysis, showed as a function of the effective frequency ( √ ν i · ν j ) as well as the fitted SEDs at order 0 (MBB), order 1 and order 2. The bottom panel quantifies the relative difference between the SIM3 cross-spectra and the same 3 fits, as defined in Eq. (16). 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 ∆D 0 (143 × 143) = 10.9 % residual, the 1 st order fit ∆D 0 (143 × 143) = 3.2 %, the 2 nd order fit ∆D 0 (143×143) = 0.5 % and the 3 rd order fit ∆D 0 (143×143) = 0.06 %.
Let us now consider a future CMB B-mode experiment that has the same frequency channels and a SNR for the dust Bmodes similar to those of Planck-HFI for the dust intensity (e.g. LiteBIRD, Hazumi et al. 2019). We know from Planck Collaboration et al. (2016b) that the dust B-mode 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 r-equivalent amplitude at 150 GHz r D (Planck Collaboration et al. 2016b), it corresponds to r D = 1.8. Our results thus 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 the dust B-modes follow a MBB SED, ∆r = 0.032 · r D = 0.07 by assuming an order 1 moment expansion, ∆r = 0.005 · r D = 0.009 by assuming an order 2 moment expansion and ∆r = 0.0006 · r D = 0.001 by assuming an order 3 moment expansion (as we have seen, parametrizing a dust decorrelation amounts to fitting the 1 st -if the parametrization is correct).
The region of the sky observed by the BICEP/Keck experiment has r D = 0.11 (BICEP2 Collaboration et al. 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 to ∆r = 0.032 · r D = 4 × 10 −3 , a second order analysis to ∆r = 0.005 · r D = 5 × 10 −4 and a third order analysis to ∆r = 0.005 · r D = 6 × 10 −5 . Mean SED in MJy 2 ·sr −1 for the 100 SIM3 crossspectra, as a function of the effective frequency ( √ ν i · ν j , in GHz), for the multipole bin centered at = 80. The MBB (solid black), order 1 (dashed yellow), order 2 (dashed-dotted blue) and order 3 (dotted red) best fits are displayed and not distinguishable. Bottom panel. Relative difference in % between the SIM3 SED and the MBB (black diamonds), order 1 (yellow diamonds), order 2 (blue diamonds) and order 3 (red diamonds) best fits.
Although these ∆r values are rough estimates (that might be overestimated in the BICEP/Keck case, because fewer SED spatial variations could occur on this small region), they give an insight on the order of magnitude of the potential bias. Therefore, they strongly advocate for the need to take into account the spectral departures from the dust MBB in future CMB B-modes analyses, targeting r values down to 10 −3 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 M 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 −3 , we need to consider all moments with an absolute amplitude larger than 10 −3 . As it can be seen in Fig. 4, most of the moments in the PR3 data set decomposition up to order 3 have an absolute amplitude larger than this threshold. In that sense, dust-dominated angular scales of the intensity PR3 data set additionally stress the need for an order 3 expansion of the polarized dust SED, in order to reach the accuracy targeted by future CMB B-mode experiments.

Summary and conclusions
In this paper, we have presented 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 of the dust SED on the sky and along the line of sight, to provide an astrophysically motivated description at the power spectrum level. The model formalism relies on the 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; they inevitably appear due to averaging effects along the line-of-sight, within the beam, and, most relevant here, due to 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 follow.
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 is 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 moments 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-modes 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 −3 .
If our results extend to polarized emission without any additional complexity, we anticipate that moment expansion up to order 3 would be required to model the dust polarization SED to the accuracy of future CMB B-modes 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 dust-dedicated frequency channels are needed in order to perform an order 2 moment expansion fit and five for order 3.
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 Bmodes frequency-dependent. Further complexity may arise due to the variation of the 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 a 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. Of course, when dealing with the polarization, other details should be carefully considered and added, 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 the extractions of primordial CMB spectral distortions. The expected signals are small (e.g., Chluba 2016) and heavily obscured by foregrounds. Current extraction methods mainly utilize information based on the SED shapes, neglecting spatial information (Sathyanarayana Rao et al. 2015;Desjacques et al. 2015;Abitbol et al. 2017). This limitation may be overcome using the techniques described here, and thus warrants further investigation. A&A proofs: manuscript no. Dust_moments_paper_AA_astroph

Appendix A: Mask and foreground templates
In this section we give 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 as LR42 and defined in Planck

Appendix B: Cross-spectra definition, covariance and correlation matrices
We give more details here 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.

Appendix B.1: Definitions
In order to avoid the noise auto-correlation bias and to reduce the level of correlated systematics, we compute the cross-spectra from data split maps (the 2 Planck half-mission maps, HM, namely HM1 and HM2   The idea of doing the sum of the 4 cross-spectra in Eq. (B.1) is to increase the SNR of D M ν i × M ν i i j , presumably at the cost of the statistical independence between distinct crossspectra, as we will see in the following. To preserve the statistical independence between the cross-spectra, we will not adopt the definition of Eq. (B.1) but that of Eq. (B.6), in Appendix B.3, for the reasons that will be presented in Appendix B.2.
The covariance matrix C from Eq. (11), 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 C 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): This correlation matrix is displayed for the SIM1 data set in Fig. B.4, for the multipole bin centered at = 100 (the shape of R is qualitatively the same in each multipole bin). It is significantly non-diagonal, showing large correlations between numerous cross-spectra. This highlights that the cross-spectra, the way we define them in B.1 are not all independent. In the following Subsection we will see why this happens. In the next, we will change the definition of the cros-spectra in B.1 to minimize these correlations.

Appendix B.2: Expected correlations from toy model
In a high signal-to-noise ratio 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's 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 high SNR 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 B.2, reads 9 : Let's 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 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 Planck-HFI noise levels, we find that R 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 2 different cross-spectra. This can happen using 2 different data splits (as the Planck HM maps) as we will see in the following.  By applying a reasoning equivalent to that of the example above to all the cross-spectra covariances in the correlation ma-trix R, 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).

Appendix 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: 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 dof = 5. When defining the cross-spectra as in Eq. (B.6), it becomes N dof = 11, allowing us to fit the SED moment expansion of the cross-spectra up to order 3 (see Sect. 5).

Appendix C: Truncating the dust moment expansion at difference orders
In Sect. 5 we have presented results from the order 3 dust SED moment expansion. We present in this Appendix the corrected dust amplitude D A D A D , the corrected spectral index β corr 0 ( ) and the relevant moment functions M in the case where Eq. (5) is truncated at order 1 and order 2 in our 3-step fitting procedure.
These other truncating order results are displayed in Figs. C.1, C.2, C.3, C.4, C.5 and C.6 for D A D A D , β corr 0 ( ), M ω 1 ω 1 , M A D ω 2 , M ω 1 ω 2 and M ω 2 ω 2 , respectively. SIM1, SIM2 and SIM3 data sets results are stable with the truncating order, while SIM4 and PR3 do change significantly. Nevertheless, SIM4 and PR3 results are affected in a very different way. As example, β corr 0 ( ) values increase with the truncating order for SIM4, while they evolve in the opposite way for SIM3. This behavior is however not observed for the moment functions. simulations. The lines are barely distinguishable. The relative difference of the dust amplitude spectrum D A D A D , the dust spectral index β 0 ( ), the corrected spectral index β corr 0 ( ) and the dust moments functions M ab up to the 3 rd order between the SIM4 and the 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 as, for example, the M A D ω i moment functions, but still well within the propagated error bars (due to division by small numbers).
According to these results we can therefore conclude that the synchrotron emission has a negligiable impact on the dust moment analysis for the Planck-HFI channels from 143 to 857 GHz.