Determining thermal dust emission from Planck HFI data using a sparse, parametric technique

Context: The Planck data releases have provided the community with sub-millimetre and radio observations of the full-sky at unprecedented resolutions. We make use of the Planck 353, 545 and 857 GHz maps alongside the IRAS 3000 GHz map. These maps contain information on the cosmic microwave background (CMB), cosmic infrared background (CIB), extragalactic point sources and diffuse thermal dust emission. Aims: We aim to determine the modified black body (MBB) model parameters of thermal dust emission in total intensity and produce all sky maps of pure thermal dust, having separated this Galactic component from the CMB and CIB. Methods: This separation is completed using a new, sparsity-based, parametric method which we refer to as premise. The method comprises of three main stages: 1) filtering of the raw data to reduce the effect of the CIB on the MBB fit. 2) fitting an MBB model to the filtered data across super-pixels of various sizes determined by the algorithm itself and 3) refining these super-pixel estimates into full resolution maps of the MBB parameters. Results: We present our maps of MBB temperature, spectral index and optical depth at 5 arcmin resolution and compare our estimates to those of GNILC as well as the two-step MBB fit presented by the Planck collaboration in 2013. Conclusions: By exploiting sparsity we avoid the need for smoothing, enabling us to produce the first full resolution MBB parameter maps from intensity measurements of thermal dust emission.We consider the premise parameter estimates to be competitive with the existing state-of-the-art solutions, outperforming these methods within low signal-to-noise regions as we account for the CIB without removing thermal dust emission through over-smoothing.


Introduction
The meticulous separation of the cosmic microwave background (CMB) from astrophysical foregrounds is a pivotal step in extracting cosmological information from the Planck data. For the intensity data several component separation techniques are currently in use; techniques such as SEVEM (Fernández-Cobos et al. 2012), NILC (Delabrouille et al. 2009), SILC (Rogers et al. 2016), SMICA (Delabrouille et al. 2003), and L-GMCA (Bobin et al. 2014) are specialised to recover only the CMB signal from all other emission sources. At present only Commander (Eriksen et al. 2008), a Bayesian parametric technique, can produce estimates of the CMB as well as the astrophysical foregrounds. More recently, a generalised version of NILC, known as GNILC, has been applied to Planck HFI data to recover the thermal dust emission (Planck Collaboration Int. XLVIII 2016). This is not the case for polarisation data, where an updated version of Commander, SMICA, and GNILC are all capable of producing estimates for polarised thermal dust and synchrotron emission (Planck Collaboration IV 2019). As our interest lies in the recovery of astrophysical foregrounds as well as the CMB and we begin our analysis using intensity data alone, we initially focus Parameter maps are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc. u-strasbg.fr/viz-bin/qcat?J/A+A/623/A21 on the specific problem of acquiring pure thermal dust emission maps tackled by Commander and GNILC.
Thermal dust emission is the dominant diffuse Galactic emission within Planck HFI data. Within the three highest HFI channels (353, 545, and 857 GHz) the other non-negligible contributions to the total intensity are the CMB, the cosmic infrared background (CIB), extragalactic point sources, and instrumental noise. Whilst the resolvable extragalactic point sources can be masked out, the unresolved extragalactic background emission that makes up the CIB requires a more global solution. The CIB makes two contributions to the total intensity measured at a particular frequency: a constant additive offset level and variations around this mean level known as the CIB anisotropies (CIBA). Planck Collaboration XI (2014) provide a method of determining the CIB offset level using correlations with H I data; while the CIB offset level can be subtracted, the CIBA cannot. The CIBA contribution at each frequency can be approximated as Gaussian distributed noise at Planck resolutions (Planck Collaboration XI 2014; Planck Collaboration Int. XLVIII 2016), but at scales smaller than a few arcminutes the CIBA traces structure formation at high redshifts (Béthermin et al. 2013).
Though not strictly a component separation problem, as there is only one well-defined diffuse emission source within an almost-Gaussian background, the Planck HFI data are an excellent testbed for techniques which aim to identify and characterise emission sources. In Irfan & Bobin (2018) we presented a new parametric technique that exploits the sparse nature of astrophysical emissions within the wavelet domain yet completes the model-fit within the pixel domain. Following the work of Planck Collaboration XI (2014) and Planck Collaboration Int. XLVIII (2016) we chose to use the modified black-body (MBB) thermal dust model. Fanciullo et al. (2015) present a comparison of several alternatives to the MBB model and find that none is able to fully represent the dust spectral energy distributions per unit extinction provided by Planck Collaboration Int. XXIX (2016). Therefore, without a clear successor to the MBB model across the full sky between 353 and 3000 GHz, we also chose this simple model with its relatively small number of parameters (temperature (T ), spectral index (β), and optical depth at 353 GHz (τ 353 )) required to characterise the recorded specific intensity I ν , where ν represents the frequency in GHz and B(T, ν) denotes the black-body function. It should be noted that the single spectral index MBB is not the only possible model; Meisner & Finkbeiner (2015) show that when considering frequencies both below and above 353 GHz the use of two spectral indices (one for low frequencies, one for high) yields an improvement in the fit. The latest Planck HFI data processing paper (Planck Collaboration III 2019), however, reveals an improvement in the calibration of the 100, 217, and 353 GHz maps. Therefore, it is possible that the need for a two-spectral index model was due to the calibration discrepancy between frequencies above and below 353 GHz in the previous data releases. More physically motivated dust particle models are also under consideration (Jones et al. 2017;Guillet et al. 2018), as are models that account for spectral deviations from the average MBB fit along a line of sight using Taylor expansion (Chluba et al. 2017) or, in polarisation, measurements of the 3D magnetic structure of dust clouds (Tassis & Pavlidou 2015;Poh & Dodelson 2017). The technique presented in Irfan & Bobin (2018), known as Parameter Recovery Exploiting Model Informed Sparse Estimates (premise), is flexible to advances in the field in that the algorithm can be updated to work with different models. In this paper we utilise the MBB model. The Planck Collaboration XI (2014) find that smoothing helps reduce the impact of the CIBA on the MBB fit, and so they only produce their all-sky maps of dust temperature and 353 GHz optical depth at the full 5 arcmin resolution of the Planck HFI data, whilst their spectral index map is at 30 arcmin. The Planck Collaboration XI (2014) thermal dust estimates can then be reconstructed using their MBB parameter maps. In place of global smoothing, GNILC (Planck Collaboration Int. XLVIII 2016) on the other hand, smoothes only in the regions where the CIBA dominates over thermal dust emission. As they apply their filtering in the wavelet domain, they can also choose to apply the greatest degree of smoothing at the angular scales which are most dominated by the CIBA. GNILC uses estimates for the CIB, CMB, and instrumental noise to inform the selection of which regions to smooth. After smoothing and recomposition back into the pixel domain they can fit the MBB model, or in fact any model, to determine the thermal dust model parameters.
Commander (Planck Collaboration X 2016) gives two 353 GHz thermal dust estimates, both determined from an MBB fit, but using different data sets. One of the dust estimates is produced at N side 256 with a 1 degree full width half maximum (FWHM) from a combination of the Planck LFI data, WMAP (Bennett et al. 2013), and Haslam (Haslam et al. 1982) data. The second, at N side 2048 with a 7.5 arcmin FHWM, uses only the Planck HFI data (excluding the 100 GHz observations) and only fits for the thermal dust spectral index, setting the thermal dust temperature at each pixel to the values determined from the low frequency fit. As Commander is a Bayesian method, the thermal dust fit makes use of priors: a Gaussian distribution with a mean of 1.55 and standard deviation of 0.1 for the spectral index and a Gaussian distribution with a mean of 23 K and standard deviation of 3 K for temperature. Liu et al. (2017) compare the thermal dust estimates produced by Planck Collaboration XI (2014), GNILC, and Commander (at N side 256) after all maps have been smoothed to a common resolution of two degrees. By computing the angular power spectra of the MBB parameters they identify a striking difference between the Planck Collaboration XI (2014) and GNILC spectral index power spectra at small angular scales. Although different features are visible in all three temperature and spectral index maps, Liu et al. (2017) show that the histogram distributions for the Planck Collaboration XI (2014) and GNILC parameters are far more similar to each other than to the Commander distributions.
An additional check of the MBB temperature and spectral index is their correlation. There may be a physical mechanism that produces an anti-correlation between T and β; Ysard et al. (2015) generate dust SEDs with various grain parameters, and from MBB fits to these pure dust SEDs they find that the variation of grain properties can result in an anti-correlation of MBB temperature and spectral index. However, it has also been demonstrated that the presence of noise biases the MBB fits to also produce such an anti-correlation (Shetty et al. 2009). Liu et al. (2017) reveal both negative and positive correlations across various regions of the sky for the Commander temperature and spectral indices, while the Planck Collaboration XI (2014) and GNILC temperature and spectral indices only display negative correlations. Liu et al. (2017) note that the discrepancies between Commander and both Planck Collaboration XI (2014) and GNILC are likely due to the priors imposed by Commander.
This paper is the empirical data accompaniment of Irfan & Bobin (2018); we use premise to determine the thermal dust MBB parameters from the Planck second data release (PR2), HFI data. Irfan & Bobin (2018) focuses on this same task, demonstrating the method on simulation data where the true parameter values are known. In Irfan & Bobin (2018) we also compose a GNILC-like algorithm of our own for comparison. Possessing the "oracle" parameter values allowed us to make to following observations: methods which rely on smoothing are at risk of smoothing out actual thermal dust information. Within the Galactic plane, premise and the GNILC-like algorithm performed comparably, but in low signal-to-noise regions premise could be seen to outperform the GNILC-like algorithm. Through a combination of -a less conservative CIBA filtering technique than GNILC, -averaging over super-pixels for the MBB fit, -the removal of non-sparse wavelet coefficients, premise was able to produce more accurate MBB parameters within low signal-to-noise regions than the GNILC-like algorithm. We now try this analysis on empirical Planck HFI data with the goal of producing all-sky MBB parameter maps -at full resolution; -in less time than a pixel-by-pixel MBB fit; -with competitive accuracy, compared to existing methods.
A21, page 2 of 17 In Sect. 2 we present the data used for our analysis; Sect. 3 gives an overview of the premise method; and Sect. 4 presents our results.

Data
In order to make comparisons between our method and existing works, we perform our analysis on the Planck PR2 data release.

Total intensity maps
We determine the thermal dust MBB parameters by fitting to total intensity data at 353, 545, 857, and 3000 GHz. The Planck PR2 data release provides the 353, 545, and 857 GHz observations across the full sky at N side 2048 and with a FWHM of 5 arcmin. We make use of the total intensity estimate at 3000 GHz provided by Planck Collaboration XI (2014). This map is a combination of the IRIS 3000 GHz map (Miville-Deschênes & Lagache 2005), made from IRAS data and the SDF 3000 GHz map (Schlegel et al. 1998). This combination favours the SDF representation on scales larger than 30 arcmin and the IRIS representation on scales smaller than 30 arcmin, up to the 5 arcmin resolution of the IRIS map. We use the version of the combined IRIS and SDF map on the Planck Legacy Archive 1 which has had the point sources removed and inpainted.
In order to account for the change of intensity across the instrumental bandpasses, colour corrections were calculated. For the 353, 545, and 857 GHz colour corrections of the Planck HFI bandpasses, the HFI bandpass responses were used (Planck Collaboration IX 2014). For the 3000 GHz colour corrections the IRAS/SFD band responses were used 2 . As the colour corrections depend on the thermal dust temperature and spectral index, they are calculated for each temperature and spectral index estimate within the model fitting routine.

Simulated CIB and instrumental noise
The premise method employs a filtering technique similar to that of GNILC and so requires estimates of the CIB and instrumental noise. We used the FFP8 (Planck Collaboration XII 2016) simulations for the Planck HFI instrumental noise. Gaussian noise with a median level of 0.06 MJy sr −1 (Planck Collaboration XI 2014) was used for the 3000 GHz instrumental noise map.
We also used the FFP8 estimates for the CIB at 353, 545, and 857 GHz, but we altered their mean values to match the CIB offsets measured by GNILC: CIB estimate -mean value + GNILC offset value. The GNILC CIB offset values are 0.1248, 0.3356, and 0.5561 MJy sr −1 for 353, 545, and 857 GHz, respectively (Planck Collaboration Int. XLVIII 2016). A 3000 GHz CIBA map of N side 2048 and FWHM 5 arcmin was created using the methodology detailed in Appendix C of Planck Collaboration XI (2014). We then fixed the CIB offset value of this map to the GNILC calculated value of 0.1128 MJy sr −1 in the same way as described for the HFI frequencies.

CMB
The premise method for estimating thermal dust requires the CMB to be removed from the total intensity maps. For this we use the CMB estimate of Bobin et al. (2016). 1 http://pla.esac.esa.int/pla 2 http://svo2.cab.inta-csic.es/theory/fps/

Additional data
We made use of the Green et al. (2018) map of interstellar reddening 3 to assess the accuracy of our optical depth at 353 GHz estimate (Sect. 4). Green et al. (2018) provide E(r p1 − z p1 ) measurements across 75% of the sky; E(r p1 − z p1 ) ≈ E(B − V) for R V ∼ 3.1. The map is based on stellar observations from Pan-STARRS1 (Chambers et al. 2016) and 2MASS (Skrutskie et al. 2006) data and is available as an all-sky HEALPix map at N side 2048.

Method
The motivation behind premise is the abundance of astrophysical information now available; we choose to make use of astrophysical models for the foreground components and aim to accurately determine these model parameters. The premise methodology is sparsity driven and although the data are fit in the pixel domain, the filtering and parameter refinement steps take advantage of the sparse nature of thermal dust emission within the wavelet domain. In this section we recap the premise method and detail the preprocessing steps specific to determining thermal dust emission from Planck HFI data. We encourage the interested reader to consult Irfan & Bobin (2018) for the mathematical derivation of the premise algorithm.
There are three main steps to the premise methodology: -Filtering: First we filter the data to suppress the CIBA and instrumental noise. We employ the GNILC technique of using a covariance matrix formed from CIB and instrumental noise estimates, but add to this method by exploiting the sparsity of thermal dust emission in the wavelet domain. As the covariance matrix is calculated for the N frequency by N pixel CIB plus instrumental noise matrix, any correlations across frequencies within the noise and CIB simulations are represented. -Super-pixel fit: Then we avoid the computational cost of completing a pixel-by-pixel MBB fit to the filtered data by fitting to super-pixels. The unique aspect of this step is that the super-pixel areas are selected by the premise algorithm itself to ensure that the regions where the MBB parameters vary most slowly are given the largest areas to average over. An initial super-pixel area of 128 × 128 pixels (roughly four by four degrees) is set in order to provide the first fit to the MBB model; this "first pass" provides the algorithm with a map of reduced chi-squared values, which it uses to select the final super-pixel areas. -Parameter refinement: Lastly, we determine the optimum parameter values per original resolution pixel. Using the super-pixel parameter maps as initial guesses, we evoke a gradient descent at each pixel to minimise the least-squares estimator. Once the gradient descent converges, we threshold the parameters in the wavelet domain to prevent non-sparse, noise terms from contributing to the final parameter estimates. We ensure MBB parameter maps at the original map resolution by using the empirical data, i.e. the total intensity maps containing thermal dust, CIB, and instrumental noise, as the observational data in our gradient descent. The filtered maps are only used within the super-pixel fit step. A flow diagram summarising these three main steps is presented in Fig. 1; as both the filtering and the super-pixel fit steps are required to obtain the initial MBB parameter estimates they are combined within the "Parameter Initialisation" block of the method flow diagram.

Pre-processing
We assume the CMB contribution at 3000 GHz to be negligible and so only subtract the Bobin et al. (2016) estimate of the CMB from the 353, 545, and 857 GHz data.
The raw data are in the form of HEALPix maps. As the quadtree works in two dimensions, we transform the (12×2048× 2048) pixel HEALPix vector into 12 (2048, 2048) matrices and perform the filtering and super-pixel fit steps of the algorithm 12 times, once for each area of 2048×2048 pixels. The method used to convert a HEALPix vector into twelve 2D faces is described in Appendix A.

Filtering
To suppress the impact of the CIBA on determining the MBB parameters, premise employs a modified version of the GNILC filtering technique on the CMB subtracted data. GNILC uses smoothing to suppress the CIBA within the regions of the total intensity maps at each scale within the spherical wavelet domain where the nuisance contributions (instrumental noise plus CIBA) dominate over the Galactic thermal dust signal. Although this technique relies upon the assumption that both the instrumental noise and CIBA can be approximated as Gaussian at each wavelet scale, the covariance matrix does take into account any slight correlations which occur over frequency.
For premise we apply the GNILC filtering as follows for each wavelet scale ( j): 1. We define the nuisance term as N = CIB estimate + instrumental noise estimate (see Sect. 2.2); 2. The raw data and nuisance estimates are divided into overlapping super-pixels of area 8 × 8 pixels with an overlapping ratio of 0.5. Our technique differs from GNILC in this respect as GNILC instead smoothes the data and nuisance estimates; 3. The N obs by N obs nuisance covariance matrix is calculated as 4. The covariance matrices of the binned total intensity maps are calculated for each bin as R tot = 1 bin area X × X T and then whitened: 5. The eigenvectors of the whitened R tot are calculated for each bin and ordered; 6. The Marcenko-Pastur distribution is used to select eigenvalues which deviate significantly from unity. Here our method again differs from GNILC, which uses the Akaike Information Criterion to select the eigenvalues that represent the signal subspace; 7. The selected eigenvectors (U s ) give the mixing matrix (F = R 1/2 nus U s ) which can be used to obtain the least-squares optimisation of thermal dust emission. However, we choose to add an additional factor to the optimisation performed by GNILC. We recapture thermal dust information lost through binning by adding a penalisation factor, to favour sparsity, into the least-squares minimisation. This penalisation factor takes the form of discrepancies between the raw data and binned dust signal (∆), which are sparse in the wavelet domain. The L2,1-norm is used to identify which discrepancies are sparse: where k is the pixel number; 8. It should be noted that a user-defined threshold value is required to identify significant wavelet coefficients, i.e. those high enough above the noise level to be considered signal. Our filtering uses two wavelet scales and a threshold cut-off value of 2.4.
We choose to alter the GNILC filtering as we believe that it may be too conservative: by striving to entirely remove the effects of the CIBA it may in fact be smoothing over the thermal dust signal itself. By using the Marcenko-Pastur distribution in place of the Akaike Information Criterion the premise filtering classifies larger regions of the full sky as signal-dominated. In addition, we then add a sparse penalisation factor to the optimisation resulting in filtered dust maps which fully capture thermal dust emission yet still contain a non-negligible CIBA contribution. This non-negligible CIBA remnant is prevented from contributing to the dust parameter estimation by the superpixel fit.

Large-scale CIB offset and point sources
After filtering to deal with the CIBA and instrumental noise we need to remove the large-scale CIB offset before we can fit the MBB model. To enable a direct comparison between our results and that of GNILC we simply use the CIB offsets given in Planck Collaboration Int. XLVIII (2016): 0.1248, 0.3356, 0.5561, and 0.1128 MJy sr −1 for 353, 545, 857, and 3000 GHz, respectively.
The 3000 GHz total intensity we use has already had its point sources masked and inpainted (Planck Collaboration XI 2014) so we only need to mask the 353, 545, and 857 GHz data. For this we use the 2015 Planck HFI point source masks.

Super-pixel area selection
We perform our MBB super-pixel fit on the filtered, large-scale CIB offset removed, point source masked data. It is worth noting that the GNILC and Planck Collaboration XI (2014) methods both apply the same MBB fitting procedure, which includes the calibration uncertainties of the Planck HFI maps. Although premise does not allow for calibration uncertainties within the data, we do not see a bias due to this within our results. Instead of working on the full (12 × 2048 × 2048) pixel matrices, the fit works in two dimensions, i.e. on each of the twelve 2048 × 2048 pixel areas. The procedure is as follows: A21, page 4 of 17 1. The instrumental noise + CIB covariance matrix is calculated for the full face within the wavelet domain for four wavelet scales; 2. The 2048 × 2048 face is split into patches of 128 × 128; 3. A single MBB is fit to each 128 × 128 patch; 4. The data-model residual within each patch is calculated per original pixel and transformed into the wavelet domain, again using four wavelet scales; 5. The reduced χ 2 at each of the four wavelet scales is calculated for each pixel. 6. The data are split into super-pixels according to the reduced χ 2 . If more than 10% of the pixels in a 64 × 64 patch contain reduced χ 2 values greater than 2 within the third wavelet scale, then the patch is split into four 32×32 patches. If more than 10% of the pixels in a 32 × 32 patch contain reduced χ 2 values greater than 2 within the second wavelet scale, then the patch is split into four 16 × 16 patches and so on; 7. A lower limit of 8 × 8 patches is enforced, and the algorithm is also prevented from cutting up the data into patch sizes that only contain masked data; 8. The MBB fit is then re-run on the algorithm chosen superpixels to obtain initial parameter estimates. Figure 2 shows the super-pixels chosen. The process of dividing the total intensity maps into twelve faces, filtering these faces, selecting super-pixels, and then fitting an MBB model to these super-pixels as described above is in fact repeated a total of seven times as we make use of a signalprocessing technique known as cycle-spinning. This technique, and specifically our application of it, is described in Appendix A.
The use of super-pixels introduces an interesting feature to our method caused by the distinct spatial features of the MBB parameters. While the MBB temperature and spectral index values vary gradually across the sky, the optical depth at 353 GHz is essentially the normalisation factor responsible for both the large-and small-scale features of thermal dust. This means that for a single super-pixel, for instance a 32 × 32 pixel area, the MBB temperature and spectral index between 353 and 3000 GHz may remain constant at 19 K and 1.55, respectively, whilst the optical depth traces out numerous filamentary features. Such a super-pixel fit would result in accurate estimates for temperature and spectral index, as the pixel area would be large enough to average out any residual noise and CIBA contributions, but provide a poorer estimate for optical depth as all variations within this area would be lost.
We deal with this trade-off between sensitivity to the thermal dust spectral form and an accurate depiction of the normalisation factor in the parameter refinement step.

Refinement
The refinement step of the premise algorithm is the final step and provides us with all-sky, full-resolution maps of the MBB parameters. As th e three MBB parameters suffer from degeneracies we first discuss temperature and spectral index and then attempt to break the degeneracy, with additional data, to determine the optical depth at 353 GHz.

Temperature and spectral index
The super-pixel fit produces initial estimates of the MBB temperature and spectral index. As these initial estimates are well informed we can use them to seed a gradient descent algorithm to determine the optimum temperature and spectral index per pixel. First we calculate an initial estimate for optical depth at 353 GHz by minimising the reduced χ 2 at each 5 arcmin pixel (k) between the total intensity (thermal dust plus CIB plus instrumental noise) data and our spectral model of thermal dust emission (S ). Our spectral model is formed using the quadtree estimates for temperature and spectral index produced by the super-pixel fit of the algorithm. We weight the data according to our CIBA plus instrumental noise estimates (n) where X[k] is a vector containing the total intensity measurement at pixel k at each frequency and where CC −1 are the inverse colour corrections for the four bandpasses and the factor of 10 20 is the conversion factor from power to MJy. We can also calculate a noise estimate: We then produce our initial τ 353 estimate by soft thresholding the map in the wavelet domain, setting the threshold at 2.4 times the noise estimate (δτ 353 ). We can then construct a thermal dust model at each frequency (ν) as τ 0 353 × B(T 0 , ν) × (ν/353) β 0 . The gradient descent is performed iteratively, per pixel, for both temperature and spectral index: The gradient path length (α) is computed analytically at each iteration; it is the minimum of the inverse Hessian matrix of the data-model residual (with respect to β and T ). The gradient descent is terminated upon convergence, i.e. when the two iterations differ in their parameter estimation by less than 10 −6 . After the gradient descent step we now have two all-sky maps, one of temperature and one of spectral index. These maps are at the full Planck resolutions of 5 arcmin as the gradient descent step uses the total intensity data (thermal dust plus CIB plus instrumental noise), not the filtered data from Sect. 3.2. Using the spherical wavelet transform we soft threshold these maps in the wavelet domain (leaving the coarse scale untouched) to remove the effect of any non-sparse elements in the total intensity maps. The threshold values are calculated at each of the five wavelet scales as twice the mean absolute deviation (MAD) of the temperature/spectral index. As variance cannot be calculated per pixel, we divide the data into patches of 32 × 32 pixels for the MAD calculation. Thresholding removes non-sparse information from the wavelet scales: if for a particular region and angular scale the signal is completely subdominant to the noise, then we lose the signal for that region at that resolution. We also exclude information from the first wavelet scale as this corresponds to individual pixel variations, so changes in temperature and spectral index across scales smaller than the Planck beam FWHM. These spurious variations are produced by the inclusion of the CIBA in the total intensity data used for the gradient descent.
We apply a total point source mask, made from the multiplication of the 353, 545, and 857 GHz Planck HFI masks, to the residual matrix ((X − Model)) in order to prevent the point sources from contributing to the final temperature and spectral index estimates. Therefore, temperature and spectral index values within the total point source mask remain at the initial, super-pixel estimate. The process of gradient descent followed by thresholding within the wavelet domain is then repeated to ensure a robust solution.

Optical depth at 353 GHz
The MBB parameters are heavily degenerate when fit in the presence of CIBA and noise. Both Planck Collaboration XI (2014) and Planck Collaboration Int. XLVIII (2016) try to combat this degeneracy through smoothing. The premise super-pixel MBB fit is specialised to select super-pixel sizes based on the spatial variation of the dust temperature and spectral index; therefore, the same method is less responsive to the overall normalisation factor (the optical depth).
To break the parameter degeneracies we returned to the total intensity data (thermal dust plus CIB plus instrumental noise) for additional information. We calculated our optical depth estimate at 5 arcmin using Eq. (3), but this time our spectral model is formed using the 5 arcmin estimates for temperature and spectral index produced by the refinement step of the algorithm. As before, we were able to calculate a noise estimate using Eq. (5) and to produce our final τ 353 estimate by soft thresholding this map in the wavelet domain, setting the threshold at 2.4 times the noise estimate (δτ 353 ).

Results
We aim to assess our thermal dust MBB parameters through comparisons with existing thermal dust templates and external data sets. In this section we refer to the results obtained by Planck Collaboration XI (2014) as "2013"; whilst the method dates from this time the results themselves have been reproduced here using PR2 data and it is these updated results that we use in this work. Figure 3 shows the all-sky MBB temperature and spectral index for thermal dust emission as calculated by premise, 2013, and GNILC. The differences between the parameter estimates are given in Fig. 4. Visually, for both the temperature and spectral index maps, the premise and 2013 maps look the most similar, differing noticeably from the GNILC temperature and spectral index maps at high latitudes. The globular features, clearly seen A21, page 6 of 17 in any of the difference maps which include the GNILC spectral index estimate, are discussed in detail in Sect. 4.6.

Temperature and spectral index: spatial distributions
In Fig. 5 we plot the percentage difference between the premise and GNILC temperature values against the percentage difference between the premise and 2013 temperature values; we also plot the same for the spectral indices. For Fig. 5 we calculate the percentage difference for the N side 2048 maps, but then downgrade the results to N side 256 for the plot in order to highlight trends over regions of the sky as opposed to individual anomalous pixels. A "spike" of linear correlation is seen in the temperature plot (top) in Fig. 5, which corresponds to a thin strip within the Galactic plane where the abundance of point sources prevents premise from providing accurate temperature estimates. Masking of this region removes the spike from the temperature density plot. Otherwise, there are no other regions where GNILC and 2013 consistently disagree with premise over the parameter values, implying that the discrepancies in parameter values are mainly due to the methodology used.
The mean spectral index and temperature values across the full sky are respectively 19.5 K and 1.54 for premise, 19.4 K and 1.60 for GNILC, and 19.8 K and 1.53 for 2013. To further probe the discrepancies in parameter estimates, Fig. 6 shows the distributions of the temperatures and spectral indices within various strips of Galactic latitude as calculated by premise, 2013, and GNILC. The mean parameter estimates within each latitude strip are noted on each histogram. At high Galactic latitudes (|b| > 40 • ) the GNILC temperature estimates are lower than those of 2013 or premise, whilst close to the Galactic plane the three temperature estimates are almost indistinguishable. The three spectral index estimates are also fairly close to the Galactic plane, varying by 0.06 within 10 • < b < 40 • and 0.03 within 0 • < b < 30 • . The largest differences in spectral index estimates are seen at high latitudes where the GNILC estimates favour higher values than those of 2013 or premise. Once again, the three sets of parameter estimates differ most notably within the regions where the CIBA and instrumental noise contribute most significantly to the total intensity.

Temperature and spectral index: anti-correlations
A possible check of MBB temperature and spectral index is their correlation. While it may be possible for there to be a physical mechanism which produces an anti-correlation between T and β, it is known that the presence of noise and the CIBA biases the MBB fits to also produce an anti-correlation (Shetty et al. 2009). Additionally, it is worth noting that the MBB is an approximated model of thermal dust emission. The very action of imposing an MBB fit onto a physical model not perfectly described by an MBB fit also introduces such anti-correlations. Figure 7 shows correlation coefficient maps made by dividing the full sky into patches of 64 × 64 pixels, giving just under two degrees of resolution. Each patch is colour-coded according to the correlation coefficient determined for the 4096 pixels within the area. The 2013 β estimates are presented at a resolution of 30 arcmin across the full sky, while their temperature estimates are at 5 arcmin. The premise temperature and spectral index estimates are at 5 arcmin resolution across the full sky and the GNILC temperature and spectral index estimates vary from 5 arcmin close to the Galactic plane to 22 arcmin within the lowest signal-to-noise regions. The correlation maps in the left column of Fig. 7   resolutions, whereas the maps in the right column are made from parameter maps smoothed to 30 arcmin. In the smoothed correlation maps larger patches of correlated and anti-correlated pixels appear at high latitudes (low signalto-noise regions) for GNILC than for premise, while the 2013 temperature and spectral index values can be seen to be barely correlated. Interestingly, all three plots show strong positive correlations around the Galactic centre which are unlikely to be caused by noise or the CIBA as the Galactic plane has the highest signal-to-noise ratio within the full sky. The same is seen for the right-hand maps; the main difference observed is the increase The first region ((l, b) = (60 • , 50 • )) is one of low signal-tonoise ratio. Here the 2013 T -β values occupy a larger parameter space indicating noisier measurements; as the correlation coefficient is inversely proportional to the standard deviation of the T -β values, this results in lower coefficient magnitudes. The GNILC parameter sets display the highest anti-correlation. The second region ((l, b) = (0 • , 3 • )) was selected as it is one of the high signal-to-noise regions where all three parameter sets display positive correlation. The third region ((l, b) = (200 • , −5 • )) is a typical high signal-to-noise region; all three parameter sets display negligible anti-correlation, though the GNILC set still posses the highest negative correlation coefficient. The fourth region ((l, b) = (150 • , 25 • )) was selected as all three parameter sets were seen to show strong negative correlations despite the high signal-to-noise ratio in this region. For regions like region 2 and region 4, where the signal-to-noise ratio is high yet significant positive/negative parameter correlations are seen for all three methods, it is possible that the thermal dust emission itself is less successfully categorised by the MBB model. Table 1 lists the properties of the four regions displayed in Fig. 8: the central coordinates of the region, the correlation coefficient between temperature and spectral index, and the standard A21, page 8 of 17

Optical depth at 353 GHz
The left column of Fig. 9 shows the all-sky τ 353 as calculated by premise, 2013, and GNILC. The middle column shows the difference maps between the three τ 353 estimates, and the right column shows the percentage difference maps between the three τ 353 estimates. The comparisons reveal the 2013 estimate to be noisier as, at high Galactic latitudes, the premise and GNILC τ 353 estimates are closer in value to each other than to the 2013 estimate. We also see again, from the percentage difference maps involving premise, the thin strip within the Galactic plane where the abundance of point sources prevents premise from providing accurate temperature estimates. The stellar reddening, E(B−V), due to extinction by interstellar dust provides a valuable and independent check on thermal dust column density estimates. Planck Collaboration XI (2014) make use of the SDSS quasar catalogue to determine correlations between E(B − V) and τ 353 . Since that analysis, a map of interstellar reddening (Green et al. 2018) has been made available; we make use of it and refer to it hereafter as the Green reddening map. Figure 10 presents correlation coefficient maps made by dividing the full sky into patches of 64 × 64 pixels, between τ 353 and E(B − V) for optical depth values calculated by premise, GNILC, and 2013. The grey patch within each map represents the 25% of the sky absent within the Green reddening map. All three maps show strong positive correlations within high signal-to-noise regions. We selected two circular regions for further investigation: a high signal-to-noise region centred at (l, b) = (90 • , −35 • ) displaying strong positive correlations for all three estimates, and a low signal-to-noise region centred at (l, b) = (124 • , 53 • ) displaying weak correlations. We consider radii of one and five degrees for each location. The relationships between τ 353 and E(B − V) in these regions are shown in Fig. 11. The data are in bins of 1000 data points for the five degree radius circles and 200 data points for the one degree From Fig. 11 it can be seen that at large angular scales the correlations between optical depth estimates and interstellar reddening are indistinguishably strong for all three optical depth estimates. At small angular scales, however, the signal-tonoise ratio within the region in question has a notable effect on the correlation coefficient. For high signal-to-noise regions all three τ 353 estimates show strong positive correlations, but within the low signal-to-noise regions weak positive correlations and even negative correlations can be seen. The GNILC estimate is more tightly correlated with E(B − V) on small scales in high signal-to-noise regions, but there are regions in the sky where the premise or 2013 estimates follow the reddening more closely than GNILC does. Regardless of which method is enlisted to deal with the contamination of CIBA, there are regions at low angular resolution which remain affected by this contamination. To try and distinguish between the three methods we now extend our analysis to thermal dust radiance estimates.

Radiance maps from the MBB parameters
Planck Collaboration XI (2014) have shown that thermal dust radiance, the integrated intensity over frequency, provides a cleaner way to view thermal dust estimates: the CIBA decorrelate over frequency, leaving less of a residual impact in the radiance maps than within the optical depth estimates. We calculate radiance using the gamma (Γ) and Riemann zeta (ζ) functions as in (Planck Collaboration XI 2014) where σ S is the Stefan-Boltzmann constant, h is the Planck constant, k is the Boltzmann constant, and ν 0 = 3.53 × 10 11 Hz. Figure 12 shows both the optical depth at 353 GHz and the radiances for premise, GNILC, and 2013 within the low and high signal-to-noise regions of Fig. 11. The five-by-five degree fields of view show the τ 353 estimates, while the zoomed-in one-byone degree insets show the radiance values. For the low signalto-noise region ((l, b) = (124 • , 53 • )) the GNILC τ 353 estimate can be seen to display a higher degree of smoothing and so fewer spatial features than the 2013 and premise τ 353 estimates; instead, within the high signal-to-noise region the three τ 353 estimates look far more similar. The premise τ 353 estimate has had thresholding applied to its wavelet coefficients and so appears smoother, even though features at the 5 arcmin scale are clearly visible upon inspection.
To fully interpret the radiance estimates shown in Fig. 12 we perform the same correlation coefficient analysis with E(B − V) as before on the same two regions, but this time for the radiance estimates. Figure 13 reveals, as before, that the three methods are only distinguishable from each other across small angular scales. For the low signal-to-noise region (top left of Fig. 13) there is an increase in the correlation coefficient for all three methods when compared to the τ 353 , one-by-one degree, low signal-tonoise scatter plot as the radiance calculation integrates the intensity over frequency and with increasing frequency the fractional percentage of CIBA in the total emission decreases. Figure 12 shows the 2013 radiance estimates within the low signal-to-noise region to be higher than the other two estimates; this is due to the steep scaling of the radiance with temperature (T 4+β ) and to the 2013 mean temperature estimate within this one-by-one degree region, which is 22.0 K as opposed to the 20.1/20.2 K mean temperature predicted by GNILC/ premise. For the high signal-to-noise region (bottom left panel of Fig. 13) the correlation coefficients decrease from optical depth to radiance. The signal-to-noise ratio is high enough within this region for the CIBA contamination within the optical depth estimates to no longer be the dominant factor. The τ 353 maps in Fig. 12 show that all three estimates have similar structures and resolutions for the highest values. Within this high signal-tonoise region the variations in MBB temperature and spectral index are responsible for the decreases in correlation coefficient. Figure 14 shows the temperature and spectral index estimates for each of the three methods per pixel of the one degree radius circle centred at (l, b) = (90 • , −35 • ). The GNILC temperature and spectral indices display the widest range of spectral indices and temperatures as well as a strong anti-correlation between these two parameters. This may be an example of the GNILC filtering smoothing out actual thermal dust emission as opposed to just the CIBA. The one-sigma deviation in temperature and spectral index for the 2013 estimates within the region shown in Fig. 14 are larger than those for GNILC and premise.

MBB parameters in detail
For a closer look at the differences between each of the MBB parameter estimates, Fig. 15 presents the T , β, and τ 353 maps for four well-known molecular clouds: Spider, Draco, Taurus, and Orion. Spider and Draco are lower in intensity than Taurus and Orion, and so their τ 353 maps are presented using a different colour scale. For the two high-intensity clouds it is hard to determine by eye any significant differences between the three τ 353 estimates; however, for Spider and particularly Draco the 2013 τ 353 maps appear to be slightly noisier. The GNILC β values for Draco and Spider are generally higher than the 2013 and premise estimates. For the Spider region the 2013 temperature estimates appear to be anti-correlated with the 2013 τ 353 estimate; this may be indicative of fitting an MBB in the presence of non-negligible noise. Within Draco the GNILC T and β estimates also appear to display correlations with their τ 353 map.
For the two high-intensity molecular clouds the only notable differences between the three temperature estimates are the visible effects of the wavelet thresholding applied to the premise estimates and the additional point-like features present only in the GNILC spectral index and temperature estimates. As the 2013 β estimates are smoothed to 30 armcin they show fewer features than the GNILC and premise estimates.
Overall, for the two low-intensity molecular clouds, premise and GNILC show more agreement in temperature with each other as 2013 seems noisier, while premise and 2013 show more agreement in β with each other than with GNILC. This may be due to the GNILC filtering which imposes higher degrees of smoothing for lower signal-to-noise regions; for the highintensity molecular clouds, where the GNILC resolution will be close to 5 arcmin, the premise and GNILC β estimates are far more similar to each other.
In Draco, Taurus, and Orion the GNILC β maps display some point-like features. All three methods deal with point sources differently: the 2013 methodology leaves point sources in the total flux maps while premise and GNILC mask out point sources. The premise method uses the Planck HFI point source masks and sets the MBB parameter values for pixels within these masks to the average value of their neighbouring pixels (the neighbourhood being determined by the super-pixel areas). GNILC creates its own point source masks from the Planck Catalogue of Compact Sources and inpaints the total flux maps within these masks. We display the combined Planck HFI point source mask (the 353, 545, and 857 GHz point sources masks multiplied together) in Fig. 16 to demonstrate the correlations between the location of point sources and the point-like features in the GNILC β maps. However premise also sees different spectral indices with the Draco cloud centre, though not in the form of compact clusters of low β values, so this may be indicative that the MBB model does not providing the best fit within molecular clouds. There are indications from both premise and GNILC, not only within Draco but also within the high-intensity molecular clouds of Taurus and Orion, that a spectral change accompanies or generates the absolute rise in opacity when going from the diffuse atomic to dense molecular media.

Thermal dust maps
In this section we compare the premise thermal dust intensity estimates to those produced by GNILC and 2013. These dust maps have been produced by reconstructing the MBB parameters: It should be noted that we need to multiply the MBB models by 10 20 for the conversion to MJy/sr and divide them by their own colour corrections (formed using each method's temperature and spectral index estimates). We re-apply the effects of the passbands so that these total intensity estimates may be visually compared with the observation data (Planck frequency maps), which have not been colour corrected. Figure 17 shows the difference between the premise MBB and GNILC MBB /2013 MBB intensity estimates at each frequency; as the Galactic plane has the highest thermal dust signal this region dominates these maps. At 353 GHz globules of negative differences can be seen in the premise MBB /GNILC MBB and 2013 MBB /GNILC MBB difference maps along the Galactic plane, this is also the case at 545 GHz and again, but too a lesser degree, at 857 GHz. The most prominent "globule" at all frequencies is situated just to the right of the Galactic centre. Figure 18 shows this particular region in detail within the premise 353 GHz dust estimate minus GNILC 353 GHz dust estimate, the GNILC temperature and spectral index maps, the premise temperature and spectral index maps, and the temperature and spectral index dif-ference maps. The globular feature is also present in the temperature and spectral index difference maps and is caused by GNILC discerning temperature and spectral index features around the Galactic plane that are not seen by premise. For the area within and around the Galactic plane, resolution is not an issue as both the premise and GNILC results are presented at 5 arcmin here. It is possible that these globular features along the Galactic plane relate to actual changes in the dust MBB temperature and spectral index, but as neither premise nor 2013 see such features they could also be a feature of the GNILC processing. Figure 19 shows percentage difference maps at each frequency: the absolute difference between the premise MBB and GNILC MBB /2013 MBB dust estimates divided by the GNILC MBB /2013 MBB dust map at each frequency. For all frequencies the largest differences are at high latitudes where the dust signal-to-noise is weakest and where the GNILC algorithm applies the greatest degree of smoothing. Along b = 0 the premise MBB thermal dust estimate differs from the 2013 MBB and GNILC MBB values to a larger extent than the 2013 MBB and GNILC MBB models differ. This highlights the difficulties the premise algorithm has in recovering model parameters across heavily masked regions. The GNILC globules are again visible across the Galactic plane at 353 GHz for both the premise MBB /GNILC MBB and 2013 MBB /GNILC MBB percentage difference maps.

CIB maps
Subtracting the MBB intensity estimates (calculated in the previous section) from the total intensity maps (after the G-MCA estimate of the CMB has been removed), then subtracting the large-scale CIB offsets results in maps of the CIBA plus instrumental noise, as shown in Fig. 20. These CIBA plus instrumental noise maps are smoothed to one degree FWHM. The expectation is for these maps to resemble Gaussian noise as the Poisson term of the CIBA only starts to dominate on angular scales smaller than the Planck HFI resolutions. We create a Galactic plane mask based on the premise τ 353 estimate, choosing to mask out the regions where τ 353 is more than 1.3 times greater than the median value. This mask leaves 57% of the sky still visible.
At 353 GHz both the premise MBB and the 2013 MBB residual maps are the cleanest, while the GNILC MBB residual map displays thermal dust contamination all over. At 545 GHz the premise MBB residual map shows negative intensity around the Galactic plane mask implying larger than true thermal dust estimates. At 857 GHz the 2013 MBB residual map shows the highest level of contamination, followed by GNILC MBB , and then premise MBB . The premise MBB residual map is cleaner than the GNILC MBB residual at 353 GHz and 857 GHz, and worse closer to the Galactic plane at 545 GHz.
The FFP8 CIBA plus instrumental noise simulations provide estimates against which the premise MBB , GNILC MBB , and 2013 MBB residual maps can be compared. In Fig. 21 we present the various CIBA plus instrumental noise estimates as power spectra. The power spectra are formed from the full-resolution residual maps. We use our Galactic plane mask to weight the a l,m to take account of the fact that these are partial sky power spectra. At 353 GHz, if the FFP8 simulations are taken to represent the true CIBA plus instrumental noise, then the GNILC MBB residual appears to overestimate the level of the CIBA plus instrumental noise, implying an underestimation in their thermal dust estimate at 353 GHz. On small angular scales, the FFP8 simulations, premise MBB and GNILC MBB residuals estimates of CIBA plus instrumental noise are all in agreement. On large angular scales, however, we see a consistent underprediction of CIBA plus instrumental noise for premise MBB and 2013 MBB residuals.
For an in-depth look at the CIBA plus instrumental noise estimates, we consider two of the four molecular clouds shown in Sect. 4.5: Draco and Spider. The Taurus and Orion molecular clouds are within our Galactic plane mask and so do not feature in our CIB analysis. Figure 22 shows the CIBA plus instrumental noise estimates for the three methods at frequencies 353, 545, and 857 GHz within Spider and Draco. We have applied the Planck HFI masks for each frequency as the total intensity data are used to form the CIBA plus noise estimates. Within both the Spider and Draco molecular clouds, the GNILC MBB residual has more regions of high intensity at 353 GHz than the other two residuals. For the Spider molecular cloud both the premise MBB and 2013 MBB residuals show a strong underprediction of the CIBA plus instrumental noise at 545 and 857 GHz as the outline of the Spider molecular cloud is still present. For the Draco molecular cloud no distinct features are seen; instead, the 2013 MBB residuals at 353, 545, and 857 GHz display a lower overall variance than the other two residuals.

Conclusions
We have presented the first all-sky, 5 arcmin resolution, estimates of thermal dust MBB parameters formed from the Planck PR2 353, 545, and 857 GHz maps as well as the IRAS 3000 GHz data. Our focus has been on obtaining model parameters that are as accurate as possible, and to gauge our success we have made comparisons with two other methods which produce MBB parameter maps: GNILC and the 2013 two-step pixel-by-pixel fit of Planck Collaboration XI (2014) updated to use PR2 data. The main difficulty in obtaining pure thermal dust maps is the presence of the CIBA; this is typically dealt with by using some degree of smoothing. We believe that the use of smoothing increases the risk of averaging out thermal dust signal itself along with the CIBA, and we try to avoid this technique by instead exploiting the sparse nature of astrophysical foregrounds within the wavelet domain.
The work presented in this paper is the application of our algorithm premise, as presented in Irfan & Bobin (2018), to Planck HFI data. Our method has three main steps: filtering of the total intensity maps, an MBB fit to adaptive super-pixel areas, and a refinement step to provide full-resolution parameter maps. Our filtering follows the philosophy of the GNILC filtering; we use noise plus CIBA estimates to identify noise dominated areas within the total intensity maps at each wavelet scale, and we bin our data thus performing a similar averaging-out to smoothing.  However, we recapture intensity variations at the 5 arcmin scale by adding a sparse penalisation factor to the least-squares optimisation. We also use a less stringent definition than GNILC for what constitutes a noise-dominated area. As a result we ensure that no thermal dust signal is lost within the filtering, the price of this being that a non-negligible level of CIBA makes it through our filtering process. The filtered maps are then divided up into super-pixels, the areas of which are specially chosen by the algorithm to provide the smallest pixel area within regions where the MBB parameters vary the most. The MBB fit to super-pixels helps to average out the remaining CIBA within our filtered maps and provides accurate initial estimates at the true parameter values. Lastly, we return to the total intensity data to perform a gradient descent at each full-resolution pixel, using our initial, accurate MBB parameter guesses to seed the descent. Any thermal dust signal smoothed out through the super-pixel averaging is recaptured along with, again, some low-level interference from the CIBA; therefore, we transform the maps into the wavelet domain one last time to remove the wavelet coefficients (at each wavelet scale) where the signal level is below that of the noise.
We find the mean MBB temperature and spectral index across the full sky to be 19.5 K and 1.54. Our temperature mean is in strong agreement with that of GNILC (19.4 K) and 2013 (19.8 K), while our mean spectral index is slightly lower than GNILC (1.60), but very similar to 2013 (1.53). Specific regions of large spectral index values seem to be present only within the GNILC spectral index maps and the GNILC temperature and spectral index values display the strongest degree of anti-correlation within low signal-to-noise regions.
For optical depth at 353 GHz the GNILC and premise estimates are in stronger agreement with each other than with the 2013 estimate at high Galactic latitudes. As the CIBA plus instrumental noise maps reveal that the 2013 dust estimate contain the largest CIBA contamination, the 2013 optical depth estimate is likely to just be noisier as opposed to biased. All three estimates show strong correlation with E(B − V) measurements at angular scales of 5 degrees but we cannot comment on the correlations at the one degree level as our analysis does not have the sensitivity to yield conclusive results on this angular scale. We reconstitute the MBB parameters for each method to produce thermal dust intensity estimates. It should be noted that GNILC provide their filtered dust maps for public use; these are not the estimates used in our analysis. Our interest lies in accurate representation of model parameters and so we reconstitute the MBB parameter maps for each method in order to provide additional tests on the parameters themselves. The difference maps between the dust estimates of all three methods reveal globular features within the 353 and 545 GHz GNILC maps, which we believe are due to the GNILC filtering technique as opposed to actual astrophysical features. Within low signal-tonoise regions all the estimates disagree with each other at 353, 545, and 857 GHz, though the premise and GNILC dust estimates display the largest percentage disagreement at 353 GHz at high Galactic latitudes and the strongest agreement at 545 and 857 GHz.
It is also seen that within the Galactic plane premise struggles to produce accurate parameter values due to the large number of point sources that obscure the thermal dust information. A21, page 16 of 17 As the premise parameter maps have been formed from Planck HFI data after the Planck HFI point source masks have been applied, our results should be considered in conjunction with these masks. The premise parameter maps cannot be used to extract thermal dust information within the regions masked out by the Planck HFI point source masks.
The CIBA plus instrumental noise maps appear to show thermal dust contamination at 353 GHz within the GNILC residual; this may account for the percentage difference between the premise/2013 and GNILC thermal dust estimates at 353 GHz at high latitudes. The power spectra reveal that there is a nonnegligible level of CIBA remaining in the 2013 and premise thermal dust estimates, on all angular scales for 2013 and angular scales larger than ∼10 arcmin for premise. We suggest that our MBB parameter estimates are an improvement on those of 2013 as our method removes more of the CIBA, and an improvement on those of GNILC as our method does not go as far as removing some of the thermal dust signal itself. pixel MBB fit on each 2D HEALPix face border effects are visible within the full-sky temperature and spectral index maps, not only at the level of the joins between the twelve 2D faces but also at the super-pixel boundaries. The majority of these border effects are erased by the parameter refinement step as the total flux maps do not contain such artefacts. However, within particularly low signal-to-noise regions the super-pixel border effects feature too prominently to be removed by wavelet thresholding; therefore, we invoke a method known as cycle-spinning (Coifman & Donoho 1995).
To cycle-spin data we apply a rotational shift to the data, operate a function on the shifted data, and then apply the inverse of the rotational shift thus changing the pixel location of any artefacts introduced to the solution by the function itself. If numerous spins are performed and their results averaged, then the operational artefacts can be averaged out with any effect to the true solution.
In the case of premise our initial temperature and spectral index estimates are formed from the full-sky intensity maps as follows: -the full-sky maps are broken up into 12 faces; -the faces are broken up into super-pixels; -a single MBB model is fit to each super-pix to determine T and β; -T and β estimates per face are reconstructed to form all sky maps of T and β. For cycle-spinning the above method is altered as follows: -a rotational transform is applied to the full-sky maps; -the full-sky maps are broken up into 12 faces; -the faces are broken up into super-pixels; -an MBB model is fit to each super-pix to determine T and β; -T and β estimates per face are reconstructed to form all-sky maps of T and β; -the inverse of the rotational transform is applied to the results. We applied six 2D rotational transforms to the full-sky maps, as described in Table A.1, and averaged the results of these six transforms plus the original (no transform) results together to form our initial estimates for temperature and spectral index.