Bayesian inference methodology to characterize the dust emissivity at far-infrared and submillimeter frequencies

ABSTRACT

We present a Bayesian inference method to characterize the dust emission properties using the well-known dust-|${\rm H\,{\small I}}$| correlation in the diffuse interstellar medium at Planck frequencies |$\nu \ge 217$| GHz. We use the Galactic |${\rm H\,{\small I}}$| map from the Galactic All-Sky Survey (GASS) as a template to trace the Galactic dust emission. We jointly infer the pixel-dependent dust emissivity and the zero level present in the Planck intensity maps. We use the Hamiltonian Monte Carlo technique to sample the high-dimensional parameter space (⁠|$D \sim 10^3$|⁠). We demonstrate that the methodology leads to unbiased recovery of dust emissivity per pixel and the zero level when applied to realistic Planck sky simulations over a 6300 |$\rm {deg}^2$| area around the Southern Galactic pole. As an application on data, we analyse the Planck intensity map at 353 GHz to jointly infer the pixel-dependent dust emissivity at |$N_{\rm side}=32$| resolution (1.8° pixel size) and the global offset. We find that the spatially varying dust emissivity has a mean of 0.031 MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$1\sigma$| standard deviation of 0.007 MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠. The mean dust emissivity increases monotonically with increasing mean |${\rm H\,{\small I}}$| column density. We find that the inferred global offset is consistent with the expected level of cosmic infrared background (CIB) monopole added to the Planck data at 353 GHz. This method is useful in studying the line-of-sight variations of dust spectral energy distribution in the multiphase interstellar medium.

1 INTRODUCTION

Galactic dust emission is one of the significant foregrounds in measurements of the cosmic microwave background (CMB) intensity and polarization at frequencies above |$\sim$|100 GHz (Planck 2018 results IV 2020). At this frequency range, the form of spectral energy distribution (SED) of the Galactic dust and the cosmic infrared background (CIB; Puget et al. 1996; Lagache, Puget & Dole 2005) are similar. While interesting in its own right, characterization and understanding of the Galactic dust properties are crucial for measuring the CMB B-mode polarization (Planck 2018 results XI 2020). To mitigate this foreground contribution from the measurements, it is essential to understand how dust grains contribute to the B-mode polarization. Understanding the Galactic dust emission is also critical for the reliable reconstruction of the CIB anisotropies (Planck intermediate results XLVIII 2016). The fact that the CIB and the Galactic dust share similar spectral properties makes it crucial to understand the spatial distribution and the SED of Galactic dust to separate the two emissions reliably. The CIB traces the matter distribution in the Universe and plays a crucial role in delensing the lensed B-mode contribution in the observed B-mode measurements (Larsen et al. 2016). It is also an important probe to be cross-correlated with other probes of the large-scale structure (Ade et al. 2014; van Engelen et al. 2015; Ade et al. 2016; Maniyar et al. 2019; Jego et al. 2023). Moreover, even though the CIB is weakly polarized, at the frequencies where CIB is significant, the polarization of the CIB can be a contaminant to CMB B-mode polarization (Feng & Holder 2020; Lagache et al. 2020).

Various methods have been used to study the nature of the Galactic dust emission with the Modified Black Body (MBB) as a model of the Galactic dust SED at Planck¹ frequencies. Inferring Galactic dust SED parameters at the highest Planck angular resolution with sufficient signal-to-noise ratio is also important because smoothing of resolution may reduce information associated with the Galactic dust. Planck 2015 results X (2016) use Commander method to perform Bayesian inference of the Galactic dust spectral properties at 7.5′ full-width half maximum (FWHM) angular resolution. A complementary method to do a similar task is implementing the dust-|${\rm H\,{\small I}}$| correlation. Galactic dust is correlated with |${\rm H\,{\small I}}$| 21 cm line emission of neutral hydrogen at high Galactic latitude and low-column density regions (Boulanger & Perault 1988). Hence, the Galactic |${\rm H\,{\small I}}$| emission map can be used as a tracer of the Galactic dust emission. Using Planck and IRAS data along with |${\rm H\,{\small I}}$| 21 cm observations obtained by Green Bank Telescope (GBT), Planck early results XXIV (2011) estimate the dust emissivities in 14 fields covering more than 800 |$\rm {deg}^2$| at high-Galactic latitude. Using the same formalism, Planck intermediate results XVII (2014) study the dust emissivity and its SED by fitting the MBB model over the Southern Galactic Pole (SGP, |$b \lt -25$|°). A proper estimation of the Galactic dust emission is needed to separate the CIB emission (see e.g. Planck intermediate results XLVIII 2016; Irfan et al. 2019; Lenz, Doré & Lagache 2019). Once the contribution of Galactic dust to a given frequency map is estimated, it can be subtracted to obtain the contribution of the CIB anisotropy (Planck 2013 results XXX 2014). Lenz et al. (2019) have produced the CIB maps by cross-correlating |${\rm H\,{\small I}}$|4PI data (HI4PI Collaboration et al. 2016) with the Planck intensities at 353, 545, and 857 GHz over approximately 25 per cent of the sky.

Planck does not measure the absolute sky background referred to as the zero levels of intensity, which have been fixed at the level of map-making (Planck 2013 results VIII 2014; Planck 2015 results VIII 2016). Two considerations go into determining the zero level of an HFI frequency map: the zero level of the Galactic emission and the CIB monopole. Zero level of the Galactic emission has been set using dust-|${\rm H\,{\small I}}$| correlation at high Galactic latitude where |${\rm H\,{\small I}}$| is the reliable Galactic dust tracer. The underlying assumption is that the dust emission is zero where |${\rm H\,{\small I}}$| column density is zero. The offset at 857 GHz is obtained by cross-correlating the Planck HFI data at 857 GHz with Leiden Argentine Bonn (LAB) Galactic |${\rm H\,{\small I}}$| Survey data. For other HFI frequencies, the Galactic zero level has been fixed using cross-correlation of maps at respective frequencies with the 857 GHz map (for details, see section 5.1 of Planck 2013 results VIII 2014). The estimated offsets are subtracted from the detector data at the time of map-making. After this step, the CIB monopole, estimated from the CIB model of Béthermin et al. (2012), has been added to each Planck HFI frequency intensity map.

Separation of the Galactic dust emission from the Planck frequency maps needs to take into account the proper treatment of the offset present in the maps. In the previous studies, the dust emissivities are fitted over local sky patches and variable offsets. Planck early results XXIV (2011) utilize dust-|${\rm H\,{\small I}}$| correlation to estimate emissivity properties of the dust at Low, Intermediate, and High-Velocity Clouds (LVC, IVC, and HVC) over 14-fields. The analysis in Planck intermediate results XVII (2014) estimates the emissivity at LVC and offset using dust-|${\rm H\,{\small I}}$| correlation over the patches of 15° diameter. Both works use the minimization technique. When dealing with very high dimensional parameter problems, |$\chi ^2$| minimization often results in parameter estimates that are far from the typical set (Mackay 2003). Furthermore, the efficiency of the |$\chi ^2$| minimization method is limited by the correlation between the parameters of interest and becomes inefficient for the high dimensional parameter space. Planck 2013 results XI (2014) measure global offset first and then estimate the Galactic dust SED parameters by fitting the MBB spectrum to offset-corrected intensity maps. In this work, we utilize the dust-|${\rm H\,{\small I}}$| correlation to jointly infer the dust emissivity and the global offset of the whole sky region considered instead of individual smaller sky patches. We use GASS |${\rm H\,{\small I}}$| data with Planck intensity maps over the same sky region used in the study of Planck intermediate results XVII (2014). We show that such an analysis can benefit from the joint inference of the emissivity and the global offset.

We sample the joint posterior distribution of emissivity and the global offset. The total number of variables to sample in this problem is around |$2 \times 10^{3}$|⁠. We use the Hamiltonian Monte Carlo (HMC) sampling method, which can more efficiently sample such large dimensional parameter space than the Metropolis–Hastings (MH) algorithm (Duane et al. 1987). For both HMC and MH algorithms, a new proposed point (⁠|$\theta ^{\star })$| is generated from the present point (⁠|$\theta)$| by taking some sort of steps based on the target density. For a D-dimensional problem, the number of steps required to reach a nearly independent proposal point grows as |$D^{1/4}$| for HMC and D for MH algorithm. The total amount of computation time with a reasonable acceptance rate grows as |$D^{5/4}$| for HMC and |$D^2$| for MH algorithm (Creutz 1988; Neal 2012; Betancourt & Girolami 2019). HMC uses the gradient of the posterior distribution to generate a proposed point, which, in principle, has an acceptance probability equal to one. This feature has led to HMC being increasingly employed in the high dimensional inference problems encountered in cosmology (Hajian 2007; Taylor, Ashdown & Hobson 2008; Jasche et al. 2010; Jasche & Wandelt 2013; Anderes, Wandelt & Lavaux 2015), including the inference of CMB foreground parameters (Grumitt, Jew & Dickinson 2020).

This paper is structured along the following lines. In Section 2, we present the data model, the likelihood, and the details of the HMC sampler. Section 3 describes the data used in this paper and the pre-processing of the data before the main analysis. Validation of the Bayesian inference method using simulated maps is discussed in Section 4. In Section 5, we present and discuss the CMB-subtracted Planck 353 GHz intensity map analysis results. We summarize the main results of the paper in Section 6.

2 METHOD

This section discusses the data model and likelihood analysis to sample the model parameters from their posterior distributions using the HMC sampler.

2.1 Dust emission model and the data likelihood

We are interested in characterizing the dust emission properties in the diffuse interstellar medium over the frequency range covered by the Planck HFI maps (⁠|$\nu \, \ge \, 217$| GHz). The major contributors to the total emission in this range of frequencies are CMB, Galactic dust, CIB, and instrumental noise. We assume that the CMB intensity map derived using the well-established component separation techniques (SMICA, NILC, SEVEM, and COMMANDER) from the Planck multifrequency observations is accurate enough in the sky region not obscured by the Galactic disc (Planck 2018 results IV 2020). With this assumption, we work with the CMB-subtracted Planck frequency maps to reduce the total number of components that need to be modelled. The |${\rm H\,{\small I}}$| column density map can be used as a tracer of the dust emission due to strong dust–gas correlation in the diffuse ISM (Boulanger & Perault 1988). We model the CMB-subtracted intensity map (⁠|$d_{\nu }$|⁠) at frequency |$\nu$| as the addition of the |${\rm H\,{\small I}}$|-correlated dust emission (⁠|$I^{{{\rm H\,{\small I}}\text{-d}}}$|⁠), the noise and a constant global offset. In general, the noise term can consist of the CIB emission |$(I^{\rm CIB}_{\nu })$|⁠, the instrument noise (⁠|$I^{\rm N}_{\nu }$|⁠), and the Galactic residual emission (⁠|$I^{\rm R}_{\nu }$|⁠). The Galactic residuals contain the dust emission from |$\tt {H}_2$| gas, which is uncorrelated with the |${\rm H\,{\small I}}$|-emission. The global offset (⁠|$O_{\nu }$|⁠) over the regions being analysed includes the contribution mainly from the CIB monopole and emission not accounted for by the |${\rm H\,{\small I}}$| emission.

It is evident from previous studies that the dust emissivity varies smoothly over the sky (Planck intermediate results XVII 2014). This fact allows us to assume the emissivity to be constant over a set of pixels. This set of pixels is determined using the healpix grid at a coarser resolution (i.e. lower |$N_{\rm side}$|⁠) than the resolution of the data map. The bigger pixel over which emissivity is assumed constant is termed a superpixel, whereas pixels within the superpixel are called subpixels. This model is expressed in the following equation,

where |$\Omega ^j_i$| indicates the direction of |$i^{th}$| subpixel within |$j^{th}$| superpixel. We model the |${\rm H\,{\small I}}$|-correlated dust emission |$(I^{{{\rm H\,{\small I}}\text{-d}}}_{\nu })$| as

where |$N_{\rm H\,{\small I}}$| denotes the integrated column density of |${\rm H\,{\small I}}$| component in the direction |$\Omega ^j_i$|⁠, and |$\epsilon ^j_{\nu }$| is the dust emissivity at frequency |$\nu$| at |$j^{th}$| superpixel. Here, we consider a single |${\rm H\,{\small I}}$| template to trace the |${\rm H\,{\small I}}$|-correlated dust emission. Because we are treating |$I^{\rm CIB}_{\nu }$| and |$I^{\rm R}_{\nu }$| as noise along with |$I^{\rm N}_{\nu }$|⁠, we call the remaining contribution of |$I^{{{\rm H\,{\small I}}\text{-d}}}_{\nu }$| and |$O_{\nu }$| as the signal model, |$s_{\nu }$|⁠,

We assume that the three noise components, the CIB, the Galactic residuals and the instrument noise, are independent. Hence, the covariance matrix of the total noise |$(\boldsymbol{\Sigma _{\nu }})$| at frequency |$\nu$| is the sum of the covariance of the individual noise component. We denote the angle between |$i^{\rm th}$| subpixel of |$j^{\rm th}$| superpixel and the |$i^{\prime \rm th}$| subpixel of |$j^{\prime \rm th}$| superpixel by |$\theta ^{jj^{\prime }}_{ii^{\prime }}$| and the elements of the covariance matrix are denoted by |$\Sigma ^{jj^{\prime }}_{\nu , ii^{\prime }} \equiv \Sigma _{\nu }(\theta ^{jj^{\prime }}_{ii^{\prime }})$|⁠. While we consider the correlation between the subpixels within a superpixel, we neglect the correlation between subpixels that belong to different superpixels, that is |$\Sigma ^{jj^{\prime }}_{\nu , ii^{\prime }} = \Sigma ^{jj}_{\nu , ii^{\prime }} \delta _{jj^{\prime }}$|⁠. Hence, the non-zero elements of |$\boldsymbol{\Sigma _{\nu }}$| are given by

where |$\boldsymbol{\Sigma _{\nu }^{\rm N}}$| denotes the instrument noise covariance, |$\boldsymbol{\Sigma _{\nu }^{\rm CIB}}$| and |$\boldsymbol{\Sigma _{\nu }^{\rm R}}$| denote the contribution to |$\boldsymbol{\Sigma _{\nu }}$| due to the CIB and the Galactic residuals, respectively. Unlike instrument noise, |$I_{\nu }^{\rm CIB}$| and |$I_{\nu }^{\rm R}$| are spatially correlated signals. Hence, its contribution to the total noise covariance matrix gives rise to the non-zero off-diagonal terms in |$\boldsymbol{\Sigma _{\nu }}$|⁠. We further assume the instrument noise, the CIB and the residual Galactic emission to be Gaussian with their respective covariance. Hence, the joint likelihood of all the data elements given the model parameters |$(\epsilon ^j_{\nu }, O_{\nu })$| is

where |$\chi ^2$| is

|$[\Sigma ^{-1}]^{jj}_{ii^{\prime }}$| represents the element of the inverse of the covariance matrix |$\boldsymbol{\Sigma }$|⁠.

In the model fitting, templates |$N_{\rm H\,{\small I}}\left(\Omega ^j_{i}\right)$| and the noise variance |$\Sigma _{\nu }$| are known, and |$\lbrace \epsilon ^j_{\nu }, O_{\nu } \rbrace$| are the unknown parameters of the model that we aim to infer from the observed data. Bayes theorem allows us to write the posterior probability distribution (⁠|$\mathcal {P}(\lbrace \epsilon ^j_{\nu }, O_{\nu } \rbrace | \lbrace d_{\nu }\left(\Omega ^j_{i}\right)\rbrace)$|⁠) of the parameters given data as

where |$\mathcal {P}(\lbrace \epsilon ^j_{\nu }, O_{\nu } \rbrace)$| is the prior probability distribution of the parameters, and |$\mathcal {P}(\lbrace d_{\nu }\left(\Omega ^j_{i}\right)\rbrace)$| is the evidence. We assume a uniform prior for all the parameters of interest without any bounds. Because the model is linear in parameters and we assume a Gaussian likelihood, the posterior distribution is also a Gaussian as a function of the model parameters. |$\mathcal {P}(\lbrace d_{\nu }\left(\Omega ^j_{i}\right)\rbrace)$| acts as a normalization constant. Hence, the functional dependence of the posterior distribution on parameters is the same as that of the likelihood distribution up to a proportionality constant. We sample the posterior distribution given in equation (7) to get the joint samples of all the parameters of interest |$\lbrace \epsilon ^j_{\nu }, O_{\nu } \rbrace$|⁠. We have around |$2 \times 10^3$| emissivity parameters per |$N_{\rm H\,{\small I}}$| template and one offset parameter.

2.2 Additional terms in the modelling

In this section, we discuss additional terms that may be required in modelling the data, their motivation, and implementation in the inference methodology.

2.2.1 Dipole term

We can model and fit for a dipole contribution in the signal model,

where |$D_{\nu }\left(\Omega ^j_{i}\right)$| accounts for residual dipole due to the CMB dipole, the CIB dipole, and the dipole from the Galactic residuals. In harmonic space, the expression for the dipole is

where |$Y_{1,m}$| are spherical harmonics with |$a^{\nu }_{1,m}$| being spherical harmonic coefficients. Using the relations |$Y_{l,-m} = (-1)^m Y^{*}_{l,m}$| and |$a_{l,-m} = (-1)^m a^{*}_{l,m}$|⁠, the above expression can be rephrased in terms of |$m = 0$| and |$m = 1$| coefficients as:

where superscripts R and I indicate real and imaginary parts of a complex quantity, respectively. In the inference process, it is convenient to treat the dipole in harmonic space because with the Gaussian likelihood, posterior is also Gaussian as a function of |$a_{1,m}$| due to their linear nature, unlike the real space variables indicating dipole amplitude and the direction of the dipole.

2.2.2 Multiple templates

The dust emission in far-infrared and sub-millimetre frequency bands can also be modelled as a linear combination of multiple |${\rm H\,{\small I}}$| templates with different dust emissivity per template. For example, Planck early results XXIV (2011) estimate the dust emissivity properties associated with LVC, IVC, and HVC clouds over 14 fields. Ghosh et al. (2017) and Adak et al. (2020) estimate the mean dust emissivity over Southern and Northern Galactic pole regions, respectively, by correlating the CMB-subtracted Planck 353 GHz map with |${\rm H\,{\small I}}$| column density associated with cold, lukewarm, and warm neutral medium (CNM, LNM, and WNM). Lenz et al. (2019) use the individual spectral channel map of |${\rm H\,{\small I}}$| brightness temperature from the HI4PI survey (HI4PI Collaboration: et al. 2016) to model the dust emission at Planck HFI frequencies |$\nu \ge 353$| GHz. The modelling and inference framework used in this work can be extended to include multiple |$N_{\rm H\,{\small I}}$| templates with the signal modelled as

where the summation is over |$N_t$| number of |$N_{\rm H\,{\small I}}$| templates (⁠|$N_{\rm H\,{\small I}}^{t}$|⁠), indexed by t, and |$\epsilon ^{j,t}_{\nu }$| is corresponding emissivity.

We test the impact of additional terms like dipole term or multiple |${\rm H\,{\small I}}$| templates on the inferred |$\lbrace \epsilon ^{j,t}_{\nu }, O_{\nu } \rbrace$| parameters using simulated maps at Planck frequencies in Section 4.

2.3 CIB model power spectra

The CIB is a relic emission from stellar-heated dust within galaxies formed throughout cosmic history. At far-infrared/sub-millimetre bands and the resolution of Planck, the CIB appears as a diffuse background emission in the total intensity. The CIB anisotropies are found to be correlated across the frequencies and follow approximately |$\ell ^{-1}$| power law angular power spectrum (Planck early results XVIII 2011; Planck 2013 results XXX 2014). We adopt the best-fitting CIB model power spectra (including the shot noise) at 217, 353, 545, and 857 GHz obtained by Planck 2013 results XXX (2014). Fig. 1 presents the model CIB power spectra used in our analysis.

We treat the CIB anisotropies as Gaussian and correlated noise in the Planck intensity maps. To take into account the correlation between pixels at a given smoothing scale, we compute the CIB covariance matrix, |$\Sigma _{\nu }^{\rm CIB}$|⁠, at a given frequency |$\nu$| between two pixels i and |$i^{\prime }$| using the relation

where |$C_{\ell }^{\rm CIB}$| is the CIB power spectrum, |$B_{\ell }$| is the beam window function of the Gaussian beam, |$P_{\ell }$| is the Legendre polynomial of order |$\ell$| and |$w_{\ell }$| is the healpix pixel window function.

2.4 Hamiltonian Monte Carlo sampling

In this section, we present the essentials of the Hamiltonian Monte Carlo sampling formalism to draw the samples of |$\lbrace \epsilon ^j_{\nu }\rbrace$| and |$\lbrace O_{\nu } \rbrace$| from the distribution given in equation (7), which is the same as the likelihood in equation (5) for a uniform prior on parameters. Details of the HMC sampling method and some considerations that went into analysing the parameter samples are discussed in Appendix  A.

HMC uses Hamiltonian dynamics to generate the proposed point and traverse the parameter space. The method treats the parameters |$\lbrace \epsilon ^j_{\nu } \rbrace$| and |$\lbrace O_{\nu } \rbrace$| as position variables and augments them with the momentum variables |$(p^{j,t}_{k\nu },\, p_{O_{\nu }})$| to define phase space dynamics.The Hamiltonian of the dynamics for the given problem is

where |$\mathcal {P}(\lbrace \epsilon ^{j,t}_{\nu }, O_{\nu }\rbrace)$| is the parameter posterior as obtained in equation (7) and |$\mu ^{j,t}_{\nu }$|⁠, |$\mu _{O_{\nu }}$| are mass terms for |$\epsilon ^{j,t}_{\nu }$| and |$O_{\nu }$|⁠, respectively. Up to a constant, which is independent of parameters of interest, |$\ln [\mathcal {P}(\lbrace \epsilon ^{j,t}_{\nu }, O_{\nu }\rbrace)]$| is

where |$\chi ^2$| is given by equation (6).

While considering the total covariance matrix, only considering the diagonal part leads to underestimating the parameter uncertainty. Whereas considering all the elements, including the ones corresponding to two subpixels of different superpixels, drastically increases the computation cost. We take the approach somewhere in between. We consider the correlations between subpixels corresponding to a given superpixel only and neglect the correlations among the pixels of two different superpixels. This is expressed in equation (4). The mathematical expressions given in the subsequent discussion are under this approximation.

We need the time derivatives of position and momentum variables to simulate the Hamiltonian dynamics. The time derivative of the position corresponding to a given parameter is simply momentum divided by the mass of the corresponding parameter:

which is not at all a computationally involved operation. The time derivative of the momentum involves computing the derivative of the logarithm of the posterior:

Our model is linear in the parameters of interest, and the likelihood is a Gaussian distribution. Hence, with the flat priors, the derivatives of the logarithm of the posterior are simple expressions given below. The derivative with respect to emissivity is

and derivative with respect to offset |$O_{\nu }$| is

While sampling the model parameters, including the residual CMB dipole contribution, the spherical harmonic coefficients |$a^{(R/I), \nu }_{1m}$| in equation (10) are jointly sampled along with emissivity and offset as |$\lbrace \epsilon ^{j,t}_{\nu }, \, O_{\nu },\, a^{(R/I), \nu }_{1m}\rbrace$|⁠. Since |$a_{1,m}^{\nu }$| in real space indicates the dipole amplitude and direction, they are sampled as global parameters similar to |$O_{\nu }$|⁠. The momentum derivative for the dipole coefficients is given by

where |$(\pm)$| corresponds to the real (R) or imaginary |$(I)$| parts respectively, of the coefficients (see equation 10). In general, |$s_{\nu }(\Omega ^j_{i^{\prime }})$|⁠, is given by equation (11). Our HMC formalism takes into account pixel-dependent dust emissivity, global offset, and three dipole amplitudes.

As the algorithm requires, we simulate the Hamiltonian dynamics using the Leap-Frog scheme (see Appendix  A). In the Leap-Frog scheme, the step size |$\Delta$| decides the time-step, which is generally different for different parameters. A general practice is to standardize the parameter distribution scales, which requires some knowledge of the parameter covariance structure (e.g. Betancourt et al. 2017). In the particular case of Gaussian posterior and the model that is linear in parameters, one can choose the mass matrix to achieve this goal. For problems where the curvature is isotropic and constant, such as for the Gaussian likelihood we consider in this work, a parameter independent |$\Delta$| can be chosen. This choice of |$\Delta$| in the case of Gaussian likelihoods is discussed in Taylor et al. (2008). For the general distributions with hierarchical modelling or non-linear parameter dependencies, this procedure may not work as well as it does in our case. By setting the mass matrix equal to the inverse of the covariance matrix of the parameters, |$\Delta$| is made independent of the distribution of the individual parameter. The inverse of the parameter covariance matrix is the negative of the parameter Fisher matrix. With our choice of neglecting the correlations between the superpixels, for the given likelihood, the Fisher matrix turns out to be diagonal over the |$\epsilon ^{j,t}_{\nu }$| parameters. The only non-zero off-diagonal terms are those which connect |$\epsilon ^{j,t}_{\nu }$| with |$O_{\nu }$|⁠, |$a_{1m}$|⁠, and |$O_{\nu }$| with |$a_{1m}$|⁠. However, we neglect these off-diagonal terms. Considering only the diagonal elements while assigning mass for the parameters then does not lead to |$\Delta$| for |$O_{\nu }$| and |$a_{1m}$| being the same as that of |$\lbrace \epsilon ^{j,t}_{\nu }\rbrace$|⁠. Hence, we choose a different |$\Delta$| in the dynamical equations corresponding to |$O_{\nu }$| and the dipole parameters. Hence, we have different step sizes in the Leap-Frog scheme, |$\Delta _{\epsilon }$| and |$\Delta _{O}$| corresponding to |$\epsilon ^{j,t}_{\nu }$| and |$O_{\nu }$|⁠, respectively. The details of this choice are discussed in Appendix  A. This leads to the mass matrix with the following terms in the diagonal. The mass for |$\epsilon ^{j,t}_{\nu }$| is

for |$O_{\nu }$|⁠, it is given by

and the following is the mass for the dipole coefficients

Note that, the mass matrix elements for |$\epsilon ^{j,t}_{\nu }$| depend on the noise covariance as well as the templates, whereas those for |$O_{\nu }$| and |$a_{1m}$| depend only on the noise covariance. In HMC, the proposed parameter is obtained by evolving Hamilton’s equations in a certain number of Leap-Frog jumps N. The product of |$\Delta$| and N determines the total distance traversed in the parameter space and controls the correlation length in the parameter chain. In this section, we have discussed the choices for N and |$\Delta$| to simulate the HMC process. While we discussed these choices without much rigour, we tested that the algorithm works using realistic simulations of the data. We would like to point out that formal methods have been developed to tune the HMC algorithm to facilitate appropriate choices for step size and path length. For example, Hoffman & Gelman (2014) present the No-U-Turn Sampler scheme to alleviate the need for the user to choose the number of steps and also presents a method for adaptive stepsize. Recent developments in also include SNAPER-HMC for implementation on GPU and TPU hardware (Sountsov & Hoffman 2021), ChEES-HMC (Hoffman, Radul & Sountsov 2021), and various adaptive schemes, for example, MALT-HMC (Riou-Durand et al. 2023). Some of these schemes are implemented in probabilistic programming frameworks such as STAN (Carpenter et al. 2017), PyMC (Salvatier, Wiecki & Fonnesbeck 2016), and pyro (Bingham et al. 2019).

We first validate the above methodology on the Planck simulations. The results obtained from the simulations are presented in Section 4. The next section discusses the data, the CIB model, and the sky masks used for our analysis.

3 DATA SETS

In this section, we describe the Planck data, |${\rm H\,{\small I}}$| data, and the sky mask used in the analysis. We also describe the procedure for computing the CIB covariance matrix using the model CIB power spectrum.

3.1 Planck data

We use the publicly available Planck 2018 Public Release 3 (PR3²) legacy intensity map at 353 GHz (Planck 2018 results I 2020) for our analysis. The 353 GHz intensity map has been provided in healpix³ (Górski et al. 2005) grid at |$N_{\rm side}=2048$| (pixel size 1.7|$^{\prime }$|⁠) with an angular resolution of FWHM 4.82|$^{\prime }$| (Planck 2018 results III 2020). We subtract the CMB contribution at 353 GHz using the following procedure. We use the |$\tt SMICA$| CMB map provided at a beam resolution of 5|$^{\prime }$| (FWHM) and |$N_{\rm side}= 2048$| (Planck 2018 results IV 2020). We smooth the 353 GHz intensity map at the resolution of the |$\tt SMICA$| CMB map using the Gaussian approximation of the Planck beam and subtract the contribution of CMB. We further smooth the CMB-subtracted Planck 353 GHz map at the beam resolution of 16.2|$^{\prime }$| (FWHM), downgrade to |$N_{\rm side}= 512$| (pixel size 6.8|$^{\prime }$|⁠) resolution. We choose the Planck 353 GHz map for our analysis because the contribution of synchrotron and free–free emissions are negligible compared to that of dust after the contribution due to the CMB is subtracted. We consider the Planck CMB-subtracted intensity map as the primary data. We treat the CIB monopole term as the global offset parameter. We use the unit conversion factors mentioned in Planck 2013 results IX (2014) to convert 353 GHz map from |${\rm K}_{\rm CMB}$| to  kJy sr|$^{-1}$| unit.

We use 300 end-to-end (E2E) noise realizations for the Planck frequency maps from the Planck Legacy Archive (Planck 2018 results XI 2020). The original maps are provided at |$N_{\rm side}= 2048$|⁠. We smoothed all the 300 noise maps at 16.2|$^{\prime }$| FWHM beam resolution and re-project them at the resolution of |$N_{\rm side}=512$|⁠. Then, we compute the variance map using the 300 smoothed E2E noise maps.

3.2 H i data

To separate the |${\rm H\,{\small I}}$|-correlated dust emission from the Planck data at high Galactic latitude, we use the |${\rm H\,{\small I}}$| data provided by GASS⁴ survey carried out by Parkes telescope (McClure-Griffiths et al. 2009). The survey observed 21 cm emission over the southern galactic sky (at declinations |$\delta \lt 1$|°) within velocity range, |$-400$| km s|$^{-1}\lt V_{\rm LSR} \lt 500$| km s|$^{-1}$|⁠; where |$V_{\rm LSR}$| is the velocity of the |${\rm H\,{\small I}}$| clouds with respect to local standard of rest. The survey has a beam resolution of 14.5|$^{\prime }$| FWHM, velocity resolution |$\delta v = 1$| km s|$^{-1}$|⁠, and root-mean-square brightness temperature uncertainty of 50 mK (⁠|$1\sigma$|⁠). The GASS survey maps used in our analysis are from Kalberla et al. (2010) and are corrected for instrumental effects, stray radiation, and radio-frequency interference.

In the southern Galactic cap, the |${\rm H\,{\small I}}$| in the Galactic disc (or what we call Galactic |${\rm H\,{\small I}}$|⁠) is mixed with significant emission from the Magellanic Stream (MS) (Nidever, Majewski & Butler Burton 2008; D’Onghia & Fox 2016). The spectra in the 3D data cube (longitude, latitude and radial velocity) likely to be associated with MS are distinguished from Galactic |${\rm H\,{\small I}}$| spectra using the velocity information (Venzmer, Kerp & Kalberla 2012). The three-dimensional model of Kalberla & Dedes (2008) helps to distinguish the spectra associated with the Galactic |${\rm H\,{\small I}}$| emission and MS. Finally, the |${\rm H\,{\small I}}$| template map is produced integrating 3D spectra over velocity range and projected on healpix grid at |$N_{\rm side}=1024$|⁠. The Galactic |${\rm H\,{\small I}}$| column density map |$N_{\rm H\,{\small I}}$| used in our analysis has an angular resolution of 16.2′ FWHM and is projected on the healpix grid |$N_{\rm side}=512$| (pixel size 6.9′). We use the Galactic |${\rm H\,{\small I}}$| map as a tracer for the |${\rm H\,{\small I}}$|-correlated dust emission. The same Galactic |${\rm H\,{\small I}}$| map is used by the Planck collaboration to study dust emission properties in the diffuse interstellar medium (Planck intermediate results XVII 2014).

3.3 Sky masking

We use the same mask as used in Planck intermediate results XVII (2014) for our dust-|${\rm H\,{\small I}}$| correlation analysis. The total sky area of the mask is 6300 |$\rm {deg}^2$| (⁠|$f_{\rm sky} = 15.3 {{\ \rm per\ cent}}$|⁠) where |$N_{\rm H\,{\small I}}\lt 6 \times 10^{20}$| cm|$^{-2}$|⁠, and thus avoids the high column density regions. The unmasked sky region covers the Galactic latitude |$b \le -25$|°. The area of 20° diameter centred around (⁠|$l_{MS}$|⁠, |$b_{MS}$|⁠)  = (⁠|$-50$|°, 0°) is masked to avoid Magellanic Stream (Nidever et al. 2010). Further, the bright radio sources at microwave frequencies and infrared galaxies at 100 μm have been masked out.

To test the dependence of the analysis results on the sky fraction, we generate four additional masks with different |$N_{\rm H\,{\small I}}$| cutoff over the range between |$2\times 10^{20} \text{cm}^{-2}$| and |$5 \times 10^{20} \text{cm}^{-2}$|⁠. We mask the regions with |$N_{\rm H\,{\small I}}$| values higher than the cutoff value. We label the mask with |$N_{\rm H\,{\small I}}$| cutoff |$Q \times 10^{20} \text{cm}^{-2}$| as MaskHIQ. The sky masks are overlapping by construction. We use overlapping masks to study the dependence of the global offset as a function of mean |${\rm H\,{\small I}}$| column density. The respective sky fraction |$f_{\rm sky}$| for each sky mask is quoted in Table 1.

4 PLANCK SIMULATIONS

In this section, we validate our methodology on the Planck simulations to simultaneously fit a dust emissivity per superpixel and a global offset using the HMC sampler.

We simulate the dust intensity maps at Planck HFI frequency bands between 217 and 857 GHz. We analyse simulated maps considering the total noise contribution from the instrument noise and the CIB anisotropies. We ignore the contribution of Galactic residuals in this analysis. We consider the instrument noise to be uncorrelated between pixels. However, for the CIB, we consider the inter-pixel correlations. Analysis with the simulated data helps to validate the pipeline and to determine some analysis choices, for example, the optimal values of the Leap-Frog step size |$(\Delta)$| and the Leap-Frog jumps (N) by examining the behaviour of the Markov Chains (MC) drawn from the posterior distribution.

We require sufficient unmasked subpixels within each superpixel to fit the signal model with the data. We set this threshold to one-third of the subpixels within each superpixel.⁵ We excluded those superpixels from the joint fitting which do not fulfil this criterion.

4.1 Simulated maps

We construct the simulated maps using the following procedure:

• We start with the dust emissivity map of Planck intermediate results XVII (2014) at 353 GHz (⁠|$\epsilon _{353}$|⁠) projected on |$N_{\rm side}=32$| (low-resolution) healpix grid. This map is obtained through the dust-|${\rm H\,{\small I}}$| correlation analysis over 15° circular patches in diameter centred on healpix pixels at |$N_{\rm side}=32$|⁠. We assume that |$\epsilon _{353}$| is the same at all subpixels defined at resolution |$N_{\rm side}=512$| that fall within a superpixel defined by |$N_{\rm side}=32$|⁠. Each superpixel contains 256 subpixels.

• We translate the dust emissivity map from 353 GHz to other HFI frequencies using the MBB spectrum,

  $$\begin{eqnarray} \epsilon _{\nu } \left(\Omega ^j_{i}\right) = \epsilon _{353} \left(\Omega ^j_{i}\right)\left(\frac{\nu }{353}\right)^{\beta _{\rm d}}\frac{B_{\nu }({T_{\rm d}})}{B_{353}({T_{\rm d}})}, \end{eqnarray}$$

  (23)

  where |$\beta _{\rm d}$| is the dust spectral index fixed to 1.5, |$B_{\nu }$| is Planck blackbody function and |$T_{\rm d}$| is the dust temperature fixed to 20 K (Planck intermediate results XXII 2015).

• We simulate the dust intensity maps at 217, 353, 545, and 857 GHz at |$N_{\rm side}= 512$| using the dust emissivity map and Galactic |${\rm H\,{\small I}}$| template. We add the global offset values from Planck intermediate results XLVIII (2016) to produce the signal model maps at all HFI frequencies.

• To simulate the instrument noise contribution, we use the variance map (⁠|$II$|⁠) computed from 300 smoothed E2E noise maps. We assume the instrumental noise to be Gaussian, white, and uncorrelated between the pixels. For the CIB noise component, we simulate the CIB map smoothed at 16.2′ FWHM beam resolution projected at |$N_{\rm side}=512$| from the model CIB power spectrum at each HFI frequency. Though the CIB is correlated between two frequencies, analysing individual frequency maps entails neglecting the correlation between the CIB emission at different frequencies.

• Finally, we co-add simulated dust intensity, global offset, instrument noise, and the CIB, all expressed in  kJy sr|$^{-1}$| units.

• To build the likelihood, we construct the instrumental noise covariance matrix (⁠|$\Sigma ^{N}_{\nu }$|⁠) and the CIB covariance matrix (⁠|$\Sigma ^{\rm CIB}_{\nu }$|⁠). |$\Sigma ^{N}_{\nu }$| is taken as a diagonal covariance matrix because we neglect the inter-pixel instrument noise. |$I^{\rm CIB}_{\nu }$| exhibits a significant correlation between subpixels within a given superpixel.

4.2 Validation with simulations

We focus our discussion of the simulated map analysis on the 353 GHz frequency channel without the loss of generality. However, we summarize the simulation results at all four HFI frequencies (⁠|$217\!-\!857$| GHz).

The output of our HMC algorithm is the chains of MC samples of the dust emissivity per superpixel at |$N_{\rm side}=32$| and the global offset. The total number of parameters at each frequency is 2011 (dust emissivity values)  + 1 (global offset) over MaskHI6.

We obtain |$2 \times 10^{4}$| samples for each derived parameter and discard the first |$10^{3}$| samples. We use the remaining |$1.9 \times 10^{4}$| MC samples for further analysis. In practice, the samples for a given parameter are not independent but are correlated with a certain correlation length. Fig. 2 presents the autocorrelation coefficients, |$\rho (t)$|⁠, for sample chains of the dust emissivity at one superpixel and the global offset chains for 353 GHz maps. The results are depicted for two choices of |$\lbrace \Delta _{\epsilon }, \Delta _{O} \rbrace$|⁠. When step size for |$\epsilon$| and |${O}$| are the same, |$\lbrace \Delta _{\epsilon }, \Delta _{O} \rbrace = \lbrace 0.1, 0.1\rbrace$|⁠, the global offset (magenta solid) chain is highly correlated. For this choice, the correlation length of emissivity and offset chain is around 8 and 220, respectively. For |$\lbrace \Delta _{\epsilon }, \Delta _{O} \rbrace = \lbrace 0.1, 1.0\rbrace$|⁠, the global offset (red dashed) chain becomes less correlated, and correlation length decreases by a factor of 2. However, the correlation length for the emissivity chain remains almost the same. For the latter choice of step sizes, the correlation length of dust emissivity and global offset chains are around 10 and 102, respectively. Therefore, we adopt the second choice of step sizes for the HMC sampling at all frequencies. Along with these |$\lbrace \Delta _{\epsilon }, \Delta _{O}\rbrace$|⁠, we choose the number of Leap-Frog steps |$N = 10$|⁠, which gives a reasonable acceptance rate and lower correlation length. We check the convergence of the chains using the Gelman–Rubin test (Gelman & Rubin 1992) (for details, see Section A1). Using seven independent chains, we find Gelman–Rubin Markov Chain Monte Carlo (MCMC) convergence diagnostics are 1.0002 (for emissivity) and 1.002 (for offset). These values confirm the chains are converged to a reasonable accuracy.

To check the correlations between the model parameters, we show the joint distributions of the dust emissivity at three superpixels and the global offset in Fig. 3 at 353 GHz. We do not find a significant correlation between the dust emissivities at two different superpixels or between the emissivity at a given superpixel and the offset. We also show the marginalized posterior distributions of the respective parameters along with their posterior mean values (red dashed line) and corresponding input values (orange dashed line) used in the Planck simulation. We find the inferred posterior mean values of the parameters agree with their respective input values.

To quantify the goodness of fit, we use posterior predictive checks considering the pixels included in the analysis, which are a subset of those in the MaskHIQ mask. This is due to the exclusion of superpixels that do not satisfy the threshold criteria for the number of subpixels. Most of the boundary pixels in the MaskHIQ mask do not satisfy the criteria and are excluded from further analysis. We term the resultant extended mask as eMaskHIQ. We replicate the simulated data based on the mean and standard deviation of the parameter samples. The distribution of the original simulated data and the replicated data match very well. Both the distributions are non-Gaussian due to the non-Gaussian nature of the thermal dust emission, which has a dominant contribution. We use summary statistics like median and 95 per cent quantile level to test the consistency between the original and the replicated simulated data. The probability-to-exceed (PTE) statistics are used to test the probability that the summary statistics of the replicated data exceed the original data. The PTE values of the summary statistics are quoted in Table 2.

From the MC chains of the parameters, we compute the mean and standard deviation of the respective parameters. In order to quantify the accuracy of inference, in Fig. 4, we present the distribution of the standard deviation (⁠|$\sigma _{\epsilon ^j_{out}}$|⁠) normalized difference, |$\delta ^j = (\epsilon ^j_{inp}-\epsilon ^j_{out})/\sigma _{\epsilon ^j_{out}}$| of input dust emissivity at all superpixels (⁠|$\epsilon ^j_{inp}$|⁠) and the posterior mean emissivity at respective superpixels (⁠|$\epsilon ^j_{out}$|⁠), for all four frequency bands. We also show the same quantity for the offset with vertical dashed lines. The mean values obtained using the MC samples agree very well with their respective mean values and are within |$3\sigma$| deviations. The largely symmetric nature of the histograms implies that the best-fitting values of output dust emissivities for all superpixels are unbiased.

In Fig. 5, we show the mean and standard deviation of dust emissivity per superpixel obtained from the dust-|${\rm H\,{\small I}}$| correlation analysis as a function of mean Galactic |${\rm H\,{\small I}}$| column density at those superpixels. The uncertainty decreases with an increase in mean Galactic |${\rm H\,{\small I}}$| column density, indicating a better estimation of dust emissivity in high column density regions within the sky mask MaskHI6. This is expected as the uncertainty scales with |$1/\mu _{O_{\nu }}$| (see equation 21) and is inversely proportional to |$N_{\rm H\,{\small I}}$|⁠.

We use the mean of the parameter samples as the best-fitting value of the respective parameter (Mackay 2003). Owing to the linear nature of the signal model, the signal that corresponds to the best-fitting values of emissivity and the offset is also the best-fitting signal. We obtain the best-fitting intensity map corresponding to the signal model using the mean values of the dust emissivity and the global offset following equation (3). In Fig. 6, we show the simulated map at 353 GHz along with maps of the best-fitting model and the residual. The residual map is the difference between the input and best-fitting model intensity maps. The residual map contains the contribution from the instrument noise and the CIB. The residual map shows the small-scale structure, lacking large-scale fluctuations, consistent with the nature of instrument noise and the CIB. While the map shows the spatial distribution of the residual, to check the nature of the distribution of the residual, in Fig. 7, we compare it (black) with the expected residual contribution from instrument noise and the CIB (magenta). We find good agreement between the recovered and the expected residual distribution at all Planck frequencies, indicating that the analysis does not introduce systematic bias. At the map level, the residual map over the unmasked region strongly correlates with the expected residual map obtained by combining the input CIB and instrumental noise map at all Planck frequencies. The Pearson correlation coefficients are 0.97, 0.96, 0.96, and 0.96 at 217, 353, 545, and 857 GHz, respectively.

Offset is a pixel-independent parameter; hence, we expect its inference to be unaffected by mask choice. We infer the global offset values for different mask choices using simulated maps to test this. In Table 1, we list the posterior mean values of the global offsets along with 1|$\sigma$| error bars and the respective masks. Irrespective of the choice of mask, the inferred offset values at all frequencies are consistent with their input values. This indicates the stability of the analysis for different choices of masks.

To elucidate the effect of the CIB noise in our analysis, we redo the analysis without the CIB contribution in the noise covariance matrix. In Fig. 8, we compare the standard deviation of the dust emissivity and offset inferred with and without considering the CIB in the noise covariance matrix. We plot the ratio of these two standard deviations for all the superpixels against the mean |$\langle N_{\rm H\,{\small I}}\rangle$| value at the respective superpixel. The uncertainty ratio in the offset parameter with and without the CIB for all four frequencies is depicted with horizontal solid lines. There is a weak dependence between the ratio of |$\sigma _{\epsilon }$| with and without the CIB as a function of |$N_{\rm H\,{\small I}}$| value. The ratio increases with the increasing frequency, consistent with the fact that the CIB noise dominates over instrumental noise at higher HFI frequency (Planck intermediate results XVII 2014).

We have shown that the method gives an unbiased inference of emissivity and the offset in the presence of realistic noise. We can faithfully take into account the noise arising due to the CIB, including the CIB-induced inter-pixel correlations within a superpixel, while we neglect the correlation between subpixels belonging to two different superpixels. The implication of considering the CIB noise is evident in Fig. 8. Not considering the CIB can lead to underestimating the uncertainty by orders of magnitude and possible bias in the inference. Further, joint sampling of emissivities with the global offset mitigates any bias that may arise due to the biased or position-dependent value of the offset.

In this simulation section, we completely ignore the contribution of the Galactic residual. Assuming the same SED for |${\rm H\,{\small I}}$|-correlated dust emission and Galactic residuals, one can expect it to dominate over the CIB at Planck’s highest frequency, 857 GHz. Like the CIB emission, we can incorporate the contribution of Galactic residuals in the noise covariance term.

4.3 Validation of two-template fit with Galactic H i and MS templates

To test the method’s robustness, we repeat the same analysis on the 353 GHz simulated map that has a contribution from the two |$N_{\rm H\,{\small I}}$| templates. We use the MS and Galactic |${\rm H\,{\small I}}$| templates to simulate map 353 GHz. The MS template traces the gas from IVC and HVC. We add the MS template with a constant emissivity |$\epsilon ^{\mathrm{MS}}_{353} = 10^{-2}\langle \epsilon ^{{\rm H\,{\small I}}}_{353}\rangle$|⁠, where |$\langle \epsilon ^{{\rm H\,{\small I}}}_{353}\rangle$| is the average over all the superpixels of the input dust emissivity map (⁠|$\epsilon _{353}$| map as used in Section 4.1). The corresponding input values for the dust emissivities are |$\langle \epsilon ^{{\rm H\,{\small I}}}_{353}\rangle = 37.3\, \,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$\epsilon ^{\mathrm{MS}}_{353} = 0.37\,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠, respectively, and the offset is |$120\,$| kJy sr|$^{-1}$|⁠, the same as that used in Section 4.1. We infer the global offset and the emissivity parameter per superpixels for both the Galactic |${\rm H\,{\small I}}$| template and the MS template over MaskHI6 performing the HMC methodology as discussed in Section 2.4. Fig. 9 shows that the recovered emissivity values are well within |$3\sigma$| of the input value without significant bias. The inferred global offset is |$120.0 \pm 0.3\, \,$|kJy sr|$^{-1}$|⁠, which is close to the input global offset value. The recovered emissivities are quoted as the mean of the dust emissivity values over all valid superpixels, along with its standard deviation. They are respectively, |$37.3\, \,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$7.2\, \,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| for the Galactic |${\rm H\,{\small I}}$| template while |$0.39\, \,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$0.27\, \,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| for the MS template.

4.4 Validation with a residual CMB dipole term

Given that we infer a global offset parameter from the partial sky, the presence of residual CMB dipole can bias its inference. We demonstrate that the method can infer a residual CMB dipole and the global offset jointly. We use this exercise to assess the extent to which the residual dipole affects the inference of the global offset from the partial sky.

We simulate a single realization of the Planck map at 217 GHz with the residual CMB dipole amplitude |$A_{\mathrm{dip}} = 3.9 \,$| kJy sr|$^{-1}$|⁠, corresponding to |$8\, \mu \mathrm{K}$| in |${\rm K}_{\rm CMB}$| unit. We choose this representative amplitude approximately equal to the difference between WMAP and Planck estimates of the Solar system dipole (Planck 2018 results I 2020). The direction |$\left(l_{\mathrm{dip}}, b_{\mathrm{dip}}\right) = \left(264.0{}^{\circ }, 48.3{}^{\circ }\right)$| is chosen same as the Solar dipole (Planck 2013 results XXVII 2014). Corresponding to these real space dipole coordinates, the input values of the harmonic coefficients in equation (10) are |$(a_{1,0}, a^{R}_{1,1}, a^{I}_{1,1}) = (5.9, 0.4, -3.7)$| kJy sr|$^{-1}$|⁠. The value of the global offset is |$O = 40\, \,$| kJy sr|$^{-1}$| (same as in Table 1), and the input dust emissivity averaged over the superpixels for 217 GHz is |$8.15\,$| kJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠. Following the procedure detailed in Section 2.4, we jointly sample emissivity per superpixel, the global offset, and the harmonic coefficients using variable Leap-Frog jumps to reduce the correlation among the parameters (as discussed in Appendix  B). We recover the offset |$O = 39.2 \pm 0.6 \,$| kJy sr|$^{-1}$| and the harmonic coefficients as |$a_{1,0}= 3.9 \pm 1.5$|⁠, |$a^{R}_{1,1} = -0.1\pm 0.3$| and |$a^{I}_{1,1} = -4.4\pm 0.4$| in kJy sr|$^{-1}$| for maskHI6. The recovered emissivity is within |$3\sigma$| as verified from a normalized histogram, indicating the method works well. The harmonic coefficients are related to the dipole amplitude and Galactic longitude and latitude, respectively, as

We calculate |$A_{\mathrm{dip}}$| and |$(l_{\mathrm{dip}}, b_{\mathrm{dip}})$| for each |$1.9\times 10^{4}$| samples using equation (24). The joint and marginalized probability distributions of the global offset, dipole amplitude, and direction are shown in Fig. 10. The corresponding mean values with standard deviations are obtained as |$A_{\mathrm{dip}} = 3.7\pm 0.3\, \,$| kJy sr|$^{-1}$| and direction |$\left(l_{\mathrm{dip}}, b_{\mathrm{dip}}\right) = \left(271.4\pm 4.2{}^{\circ }, 31.0\pm 10.9{}^{\circ }\right)$|⁠.

We find that the amplitude of the fitted dipole is positively correlated with the global offset. This is expected due to the dipole direction lying in the northern Galactic part, which is almost opposite to the SGP. The increased dipole amplitude corresponds to more negative dipole fluctuation in the SGP region, which is compensated by the global offset parameter increase. A similar reason leads to the positive correlation between the |$b_{\mathrm{dip}}$| and the offset. Lower values of |$b_{\mathrm{dip}}$| correspond to the most negative part of the dipole pointing away from the SGP, leading to the dipole compensating for a greater fraction of the positive zero level and letting the offset have a relatively smaller value. Including dipole in the analysis does not significantly change the inference of the offset, indicating unbiased inference. However, due to the correlated nature of the offset and the dipole parameters, there is a significant increase in the uncertainty of the offset parameter. Table 3 presents the recovered offset and dipole parameters as a function of sky masks at 217 GHz. As the sky area increases, the uncertainty of the global offset and three dipole parameters decreases.

We repeat the analysis with the residual CMB dipole at 353 GHz using the same Planck simulations. The same residual CMB dipole amplitude at 353 GHz is |$A_{\mathrm{dip}} = 2.4\, \,$| kJy sr|$^{-1}$| (in intensity units). We performed the same analysis and recovered the offset |$O=119.0\pm 2.4\, \,$| kJy sr|$^{-1}$| and dipole amplitude |$A_{\mathrm{dip}} =3.7 \pm 1.4\, \,$| kJy sr|$^{-1}$| and direction as |$\left(l_{\mathrm{dip}}, b_{\mathrm{dip}}\right) = \left(277.4\pm 31.1{}^{\circ }, 14.9\pm 43.2{}^{\circ }\right)$|⁠. The difference between input and the recovered values of the O, |$A_{\mathrm{dip}}$|⁠, |$l_{\mathrm{dip}}$| and |$b_{\mathrm{dip}}$| are respectively |$0.4\sigma$|⁠, |$-0.9\sigma$|⁠, |$-0.4\sigma$| and |$0.8\sigma$|⁠, |$\sigma$| being the uncertainty on the respective parameters. The uncertainty of inferred dipole parameters at 353 GHz is larger than that at 217 GHz. This is due to higher noise at 353 GHz, where the contribution from the CIB is higher than at 217 GHz (see Fig. 1). The lower amplitude of residual CMB dipole at 353 GHz (in intensity units) and higher uncertainty results in a low signal-to-noise ratio, but the recovered value is consistent with the input within |$1\sigma$|⁠. Our results show that the residual CMB dipole contribution can be ignored at frequencies 353 GHz and above. We conclude that at the noise level considered here, if the offset is larger than the amplitude of the dipole, neglecting the residual CMB dipole would not lead to significant bias in the inference of global offset.

5 PLANCK DATA RESULTS

In this section, we discuss the analysis results of the 353 GHz CMB-subtracted Planck intensity map. We model the data with pixel-dependent dust emissivity at |$N_{\rm side}=32$| resolution and a global offset. For this analysis, we do not consider the residual CMB dipole contribution in the signal model. As shown with the simulations, the noise at 353 GHz leads to increased uncertainty on the residual CMB dipole parameters as compared to the same inference at 217 GHz. However, we do test the robustness of the offset to the addition of the MS template, discussed later in this section.

We use the same superpixel and subpixel resolution as in the analysis of the simulations. We apply the HMC sampler as discussed in Section 2.4. We use the same values for the HMC sampler hyper-parameters (Leap-Frog step sizes and the number of jumps) as were used in the simulations for the analysis of the Planck data. We obtain |$2 \times 10^4$| samples of the emissivity and the offset parameters over MaskHI6, discard the first |$10^{3}$| samples, and perform the analysis with the remaining |$1.9 \times 10^{4}$| samples. Using the remaining |$1.9 \times 10^{4}$| MC samples, we compute the posterior mean and standard deviation of the model parameters. The Gelman-Rubin convergence diagnostics values for the Planck 353 GHz data estimated using seven independent chains are 1.0003 (for emissivity) and 1.003 (for offset). This shows that the chains are converged.

In Fig. 11, we show the map of the mean and the standard deviation of the dust emissivity for analysis with the MaskHI6 mask. The mean |$\epsilon$| map shows the fluctuations that extend down to the scales of superpixels, indicating small-scale variations (1.8° pixel size) in dust emissivity. The standard deviation map shows the inhomogeneity in the uncertainty of the dust emissivity. The |$N_{\rm H\,{\small I}}$| map partly determines the inhomogeneity in |$\sigma _{\epsilon }$| map as the instrument noise and the CIB are fairly statistically isotropic. This map also allows us to assess the spatial variation in the dust emissivity. The spatial average of the dust emissivity map in Fig. 11 is |$0.031 \, \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and the standard deviation of the dust emissivity values over the unmasked sky area is |$0.007\, \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠. Note that, as shown in Fig. 12, the spatial distribution of the dust emissivity is slightly skewed towards higher values. The upper panel of Fig. 12 shows the histogram of dust emissivity values in the |$\epsilon$| map depicted in Fig. 11. The lower panel of Fig. 12 shows the variation of dust emissivity as a function of |$N_{\rm H\,{\small I}}$|⁠. The data point and the error bar represent the mean and standard deviation of the emissivity computed over superpixels that fall within the given |$N_{\rm H\,{\small I}}$| bin. The average dust emissivity increases with increasing |${\rm H\,{\small I}}$| column density. Our spatial average value of the dust emissivity is |$\approx 20 {{\ \rm per\ cent}}$| lower than the value quoted in Planck intermediate results XVII (2014). This difference could arise due to (1) different aperture sizes considered for the dust-|${\rm H\,{\small I}}$| correlation analysis in Planck intermediate results XVII (2014) and (2) change in the calibration of the Planck 353 GHz map from PR1 to PR3. For MaskHI6, we report the global offset value equal to |$0.1284\pm 0.0003\, \,$| MJy sr|$^{-1}$|⁠.

We obtain the signal model map using the mean dust emissivity map and the mean global offset. As superpixels do not overlap and the emissivity varies over the superpixels, the estimated mean dust emissivity map becomes discontinuous at the scale of healpix superpixel with 1.8° angular size (pixel size at |$N_{\rm side}= 32$|⁠). To avoid this discontinuity, we first obtain the dust emissivity map at |$N_{\rm side}= 512$| by interpolating the |$N_{\rm side}= 32$| pixel values using interpolate.griddata function of scipy and then smooth the resultant mean dust emissivity map with a Gaussian kernel of 1.8° FWHM. We construct the signal model map using the smoothed dust emissivity map and the mean value of the global offset following equation (3). We show the results of this analysis in Fig. 13. The leftmost panel in Fig. 13 shows the CMB-subtracted Planck intensity map at 353 GHz over the extended unmasked region eMaskHI6, and the second panel shows the signal map. The third panel in Fig. 13 shows the difference between the CMB-subtracted Planck 353 GHz intensity map and signal model map over the same sky mask.

To compare the residual distribution with the expected distribution of the residual at 353 GHz, we simulate a random Gaussian realization of the CIB at a beam resolution of 16.2′ (FWHM) using the best-fitting CIB model angular power spectrum and an uncorrelated Gaussian white noise realization from the E2E smoothed noise variance map and add them. The expected distribution of the total noise is centred around zero because both the CIB and the instrumental noise have a zero mean. If our model assumptions are correct, the distribution of the residual should be consistent with the distribution of the expected fluctuations from the CIB and instrument noise, similar to as seen in the simulated data (see Fig. 7). The fourth panel of Fig. 13 shows the residual for masks with different sky fractions and compares it with the residual expected from the CIB and the instrument noise for the respective sky fraction. Though the residual distributions obtained using different MaskHIQ masks are mostly consistent with the expected model residuals from the CIB and instrumental noise contributions, they have a minor fraction of pixels contributing to the non-Gaussian tails with increasing sky fractions. The residual histogram is in logarithmic scale to highlight the asymmetric nature of the positive tail arising from a very small number of localized pixels in which dust is uncorrelated with the Galactic |${\rm H\,{\small I}}$|⁠.

Table 4 shows the recovered offset and recovered mean emissivities at 353 GHz over different sky masks. As the available sky area decreases, the mean value of dust emissivity decreases, indicating a |${\rm H\,{\small I}}$| column density dependence of dust emissivity. We also see an opposite trend in the global offset value, which increases with the reduction in the sky area and decreasing average emissivity. This is expected as a decrease in the global offset value compensates for an increase in the mean dust emissivity. Over the range of |$N_{\rm H\,{\small I}}$| thresholds and corresponding sky area, the inferred offset parameter ranges from 0.1358 MJy sr|$^{-1}$| (over MaskHI2) to 0.1284 MJy sr|$^{-1}$| (over MaskHI6). This range encompasses the value of the CIB monopole in the 353 GHz intensity map. The offset we infer is the total zero-level of the map that includes the CIB or any other residual zero-level. The differences between the CIB monopole value and the value we infer could be due to the differences in (1) the smoothing scale of the Planck and |${\rm H\,{\small I}}$| column density map, (2) the Planck Galactic zero-level is estimated over a larger sky fraction (28 per cent of the sky) than considered here, and (3) the difference in the local velocity cloud |${\rm H\,{\small I}}$| column density threshold (⁠|$\lt 3 \times 10^{20}$| cm|$^{-2}$|⁠) from the LAB survey (Kalberla et al. 2005).

To check the robustness of the global offset value due to the modelling error, we perform the inference with the two-template fit using Galactic |${\rm H\,{\small I}}$| and MS templates using the methodology discussed in Section 2.2.2. We infer the offset |$O = 0.1281\pm 0.0003 \,$| MJy sr|$^{-1}$| at 353 GHz over MaskHI6 for the two-template fit. The average of best-fitting values of recovered emissivities for Galactic |${\rm H\,{\small I}}$| and MS templates are |$0.031\, \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$-0.6 \times 10^{-4} \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠, respectively. The corresponding |$1\sigma$| deviation on them are |$0.007\, \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$| and |$1.8\times 10^{-4} \,$| MJy sr|$^{-1}$||$(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠, respectively.

6 SUMMARY

We demonstrate the use of the Bayesian inference method to estimate the pixel-dependent dust emissivity and the global offset in the diffuse interstellar medium at far-infrared and submillimeter frequencies. We utilize the dust-|${\rm H\,{\small I}}$| correlation to model the Galactic thermal dust emission, using the integrated Galactic |${\rm H\,{\small I}}$| column density map from the GASS survey (McClure-Griffiths et al. 2009) as a template. The signal model incorporates spatially varying dust emissivity and global offset over a 6300 deg|${}^2$| region centred around the southern Galactic pole. The HMC method allows efficient sampling of the high dimensional posterior distribution of the dust emissivity and the offset parameters.

We first validate the method on the Planck simulations, which include the CIB and instrumental noise. The dust emissivity parameters are fixed based on earlier Planck analysis (Planck intermediate results XVII 2014). This validation process allows us to test the analysis pipeline and fix the parameters of the HMC sampler that we use for the real Planck data. Given that the data are on the partial sky, the inference of the offset can be biased if any residual CMB dipole is not taken into account. We demonstrate the nature of the inferred parameters in the presence of residual CMB dipole at 217 GHz. We show that the amplitude of the residual CMB dipole is highly correlated with the global offset parameter for partial sky analysis. A small residual CMB dipole does not significantly affect the global offset inference beyond the increase in the error bar of the global offset, whose uncertainty is still significantly smaller than the global offset value at 217 GHz. At 353 GHz, the error bar on the fitted residual CMB dipole term increases due to increased noise contributions from the CIB and the instrumental noise. However, we note that the joint dipole and offset inference will be important in the case of application to the NPIPE data where the CMB dipole is retained in the frequency maps (Planck intermediate results LVII 2020). We further test the robustness of the inferred offset in the presence of multiple |${\rm H\,{\small I}}$| templates in the signal model. We consider the emission from the MS as an additional component in the signal model. The global offset value does not change significantly when considering the two-template analysis – Galactic |${\rm H\,{\small I}}$| and MS templates.

As a demonstration of the method, we apply the same HMC methodology to the Planck intensity map at 353 GHz. Results of this analysis are depicted in Figs 11, 12, 13, and Table 4. As shown in Fig. 11, we infer the emissivity of the dust over the sky area 6300 deg|$^2$| and its variation at scales larger than 1.8°. For the region of interest over 6300 deg|$^2$| area centred around the southern Galactic pole with |$N_{\rm H\,{\small I}}$| threshold |$N_{\rm H\,{\small I}}\lt 6 \times 10^{20} \text{cm}^{-2}$|⁠, the spatial average of dust emissivity is |$0.031\, \mathrm{MJy\, sr^{-1}}(10^{20} \, \mathrm{cm^{-2}})^{-1}$| with standard deviation of |$0.007\, \mathrm{MJy\, sr^{-1}}(10^{20} \, \mathrm{cm^{-2}})^{-1}$|⁠. The inferred offset for MaskHI6 is |$0.1284\pm 0.0003\, \,$| MJy sr|$^{-1}$| is close to the monopole of the Béthermin et al. (2012) CIB model value |$0.130 \pm 0.038\, \,$| MJy sr|$^{-1}$| added to the Planck intensity map at 353 GHz after setting the Galactic zero level (Planck 2018 results I 2020). We further performed the same analysis on the smaller sky mask by putting lower |$N_{\rm H\,{\small I}}$| thresholds. As expected, we find the inferred value of the global offset is stable for different sky masks. We also find that the mean value of the dust emissivity decreases as we go to the low column density regions or lower sky masks. The non-Gaussian tail in the residual distribution (caused by the emission not correlated with |$N_{\rm H\,{\small I}}$|⁠) does depend on |$N_{\rm H\,{\small I}}$| threshold and sky fraction. This additional emission component can be partly dealt with as an additional noise component in the data model, denoted by |$I^{\rm R}_{\nu }$| in equation (1), which we leave for future work.

The methodology introduced in this paper opens a new way to jointly estimate the pixel-dependent dust emissivity and a global offset over the field of interest using the dust-|${\rm H\,{\small I}}$| correlation. The HMC method is able to sample around |$2\times 10^3$| parameters and infer the parameter posterior distribution. This method can be useful in studying the 3D variation of dust SEDs to constrain the frequency decorrelation of dust B-modes. In the future study, we will apply this technique to study the dust emissivity properties associated with different |${\rm H\,{\small I}}$| phases (CNM, LNM, and WNM) of diffuse ISM considered in Adak et al. (2020) and Ghosh et al. (2017), which will be useful to extend the sky model of dust polarization at multiple sub-mm frequencies. The residual maps estimated using this method will be useful for estimating the CIB maps at the best angular resolution of Planck. In the forthcoming paper, we will apply the same methodology to other Planck HFI frequency maps to separate the thermal dust emission from the CIB emission and characterise the dust emissivity parameters.

Acknowledgement

The Planck Legacy Archive (PLA) contains all public products originating from the Planck mission, and we take the opportunity to thank ESA/Planck and the Planck Collaboration for the same. This work has used |${\rm H\,{\small I}}$| data of the GASS survey headed by the Parkes Radio Telescope, a part of the Australia Telescope funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. Some of the results in this paper have been derived using the healpix package. All the computations in this paper were run on the Aquila cluster at NISER. This work was supported by the Science and Engineering Research Board, Department of Science and Technology, Govt. of India, grant number SERB/ECR/2018/000826. DA acknowledges UGC for providing Ph.D. fellowship when a major part of this project is done. DA acknowledges IMSc for providing a postdoc fellowship for one year when the rest of the project is completed. SS acknowledges support from SERB for the Research Associate position at NISER, during which a part of this work was performed, and the Beus Center for Cosmic Foundations for current support. SS would like to thank NISER Bhubaneswar for the postdoctoral fellowship. GL acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 788212) and from the Excellence Initiative of Aix-Marseille University-A*Midex, a French ‘Investissements d’Avenir’ programme.

DATA AVAILABILITY

The Planck 353 GHz intensity map and the component-separated |$\tt SMICA$| CMB map are available at https://pla.esac.esa.int. The Galactic |${\rm H\,{\small I}}$| column density map is downloaded from https://www.astro.uni-bonn.de/hisurvey/index.php. The different sky masks, the dust model map and the residual map at 353 GHz generated in this analysis are publicly available at https://github.com/shabbir137/dust_emissivity.git.

Footnotes

References

APPENDIX A: DETAILS OF HMC IMPLEMENTATION

Different MCMC sampling methods differ from each other mainly in the way they generate the proposed point. HMC uses Hamiltonian dynamics to generate the proposed point (Duane et al. 1987; Neal 2012). This is accomplished by introducing an additional set of parameters, called momentum parameters, one momentum parameter (⁠|$p_i$|⁠) corresponding to each parameter of interest (⁠|$q_i$|⁠). The momentum variables are chosen to follow a Gaussian distribution with a covariance defined by a mass matrix. One has to choose the mass matrix specific to the problem. This aspect is similar to the choice of parameters of the proposal distribution in sampling with the Metropolis–Hastings algorithm. The multidimensional Gaussian distribution of momenta is augmented with the original probability distribution that needs to be sampled. For example, if we have ‘|$q_i$|’ as n number of parameters of interest with the probability distribution |$\mathcal {P}(\lbrace q_i\rbrace)$|⁠, the probability distribution that is sampled in HMC is the joint probability distribution of |${p_i}$| and |${q_i}$|⁠:

Here, |${\bf p}$| is the vector of momenta |$\lbrace p_i\rbrace$| associated with the parameters and |$\mathcal {H}$| is the Hamiltonian comprising of a Kinetic Energy term and a Potential Energy term:

In the above equation, the first term on the right-hand side is the kinetic energy term, and the second term, |$- \ln [\mathcal {P}(\lbrace q_i\rbrace)]$|⁠, acts as potential energy. |$\boldsymbol{\mu }$| is the mass matrix corresponding to the set of parameter |$\lbrace q_i \rbrace$|⁠. In general, |$\boldsymbol{\mu }$| is chosen as equal to the inverse of the covariance matrix of the parameters of interest. This choice may lead to |$\boldsymbol{\mu }$| having off-diagonal elements, leading to more computational cost. Instead, we choose the diagonal mass matrix. Further details about this choice and its implications are discussed later in this section. With a diagonal mass matrix, the resulting Hamiltonian is

Hamilton’s equations for variable q and the conjugate momentum p are

The above equations are solved using symplectic integration methods. In this work, we choose the Leap-Frog method to simulate the Hamiltonian dynamics.

In the above equations, |$\dot{p}(q(t))$| is the momentum derivative substituted from equations (17) or (18) depending on the parameter.

Below is the schematic of the algorithm we followed in sampling the posterior distribution. The algorithm is written in terms of the generic variables |$q_i$| (which refer to |$\lbrace \epsilon ^{j,t}_{\nu }, \, O_{\nu },\, a^{(R/I), \nu }_{1m}\rbrace$|⁠) and|$p_i$|⁠.

1: Initialize q⁽⁰⁾

2: fork = 1...N_sdo

3:  |$p \sim \mathcal {N}(0, \mu)$|⁠; where |$\mathcal {N}$| is Gaussian distribution with variance μ.

4:  |$q_{(0)}^{*}, p_{(0)}^{*} = q^{(k-1)}, p$|

5:  forj = 1...Ndo

6:   a Leap-Frog move from |$\left(q^{*}_{(j-1)}, p^{*}_{\left(j-1\right)}\right)$| to |$\left(q^{*}_{(j)}, p^{*}_{(j)}\right)$|

7:  end for

8:  |$q^{*}, p^{*} = q^{*}_{(N)}, p^{*}_{(N)}$|

9:  Accept q*, p* with the probability |$\rm {min} \left(1, e^{-H(q^{*}, p^{*}) + H(q^{(k-1)}, p^{(k-1)})}\right)$|⁠. If the proposed point is rejected, then q^k, p^k = q^{k − 1}, p^{k − 1}

10: end for

As is common in the MCMC methods, we start the algorithm with a reasonable guess of the parameter values |$\lbrace q^{(0)}_i\rbrace$|⁠. The corresponding momenta |$\lbrace p^k_i\rbrace$| are drawn from a multinormal distribution at the start of each |$k^{th}$| step, including the very first step. With this as the initial point on the phase space trajectory, Hamilton’s equations evolve these variables in the phase space. This evolution is simulated by N Leap-Frog jumps to arrive at a proposed point in the phase space |$q^{*}, p^{*}$|⁠. The leap-frog method is volume-preserving and time-reversible. However, Hamiltonian |$(\mathcal {H})$| may not be conserved in practical numerical computations due to leapfrog discretization, which can introduce a bias with respect to the target distribution. This bias can be reduced with a sufficiently small step size, which leads to less discretization error. To avoid this bias, we use the Metropolis rule, according to which the proposed point is accepted with the probability equal to the minimum of |$(1, e^{-H(q^{*}, p^{*}) + H(q^{(k-1)}, p^{(k-1)})})$|⁠.

Once the samples of |$(p_i, q_i)$| are obtained, getting the Monte Carlo samples of |$q_i$|⁠, which samples the distribution |$\mathcal {P}(q_i)$|⁠, is straightforward. Discarding the samples of |$p_i$| from the joint samples of |$(p_i, q_i)$| leads to the marginalization of the distribution |$\mathcal {P}(q_i, p_i)$| with respect to |$p_i$|⁠. As |$p_i$| and |$q_i$| are independent variables, the resultant samples of |$q_i$| are the desired samples of |$q_i$| drawn from the probability distribution of our interest |$\mathcal {P}(q_i)$|⁠.

A1 Burn-in, correlation length, and convergence test

Here, we discuss some considerations for analysing these Markov chains before using these chains to draw the inference.

A1.1 Burn-in

We neglect some of the initial samples from the chain as Burn-in. To quantitatively determine the Burn-in sample size, we monitor the |$\chi ^2$| using a model map estimated with the increasing number of samples in the chains of the parameters. The model map is built from the posterior mean of the parameter chain with an increasing number of samples. Burn-in sample then consists of the initial points where |$\chi ^2$| is away from the total number of degrees of freedom. In practice, we discard the first |$10^{3}$| samples, which is much larger than the Burn-in sample determined using |$\chi ^2$| criteria. Note that, in the absence of explicit priors on the parameters and Gaussian likelihood, |$\chi ^2$| is equal to the logarithm of the posterior up to a constant term.

A1.2 Correlation length

When we draw the samples from the target distribution, all the samples may not be independent because of the correlation within a chain of samples. The effective number of the independent samples out of the total samples, |$N_s$|⁠, is |$n_{\rm eff} = N_s/L$|⁠, where correlation length, L, is defined as,

where |$\rho (t)$| is the autocorrelation coefficient of the chain for lag t (Taylor et al. 2008).

Correlation length depends on the total jump, |$s = N \Delta$| in the Leap-Frog scheme and the scale of the distribution of the parameters of the problem. The total jump s should be of the order of scale of distribution of the parameters to get less correlated samples. |$\Delta$| controls the accuracy with which we implement Hamiltonian dynamics. This, in turn, determines the acceptance rate of the proposed sample. For a given |$\Delta$|⁠, N determines how far the proposed sample is from the current position. If |$\Delta$| is chosen to be too small, N needs to be large enough to move a sufficient distance along the trajectory, resulting in increased computation time. In contrast, for the choice of a large value of |$\Delta$|⁠, the computation of the phase space trajectories becomes relatively less accurate, resulting in a reduced acceptance rate. If N is chosen to be small, samples are correlated, whereas too large a value of N may bring back the proposed sample very close to the starting point after completing Leap-Frog jumps (Neal 2012). In Appendix  B, we show the difference in correlation length for fixed and variable Leap-Frog jumps. The accuracy of Hamiltonian dynamics can be increased by efficient choice of total jump in the Leap-Frog scheme (Hoffman & Gelman 2014; Bou-Rabee & Sanz-Serna 2017; Hoffman et al. 2021; Sountsov & Hoffman 2021).

As discussed in Section 2.4, we choose the mass matrix |$\boldsymbol{\mu }$| to make step size |$\Delta$| independent of the scale of the target distribution. In the approximations considered in our analysis, neglecting the correlation between subpixels belonging to two different superpixels renders the mass matrix sparse and diagonal dominant. The only non-zero off-diagonal terms are those due to the cross-correlation among emissivity, offset and spherical harmonic coefficients (see Section 2.4). The rest of the off-diagonal terms are the cross-correlation between emissivities at different superpixels, which are zero due to the approximations considered in our analysis. Inverting the mass matrix and evolving the vector forms of the Leap-Frog equations with |$\boldsymbol{\mu }^{-1}$| are relatively computationally costly when dealing with the non-diagonal mass matrix. Therefore, we neglect the off-diagonal terms in the column (and row) of the mass matrix that connect the emissivity and offset parameters. As a consequence of this choice, we have to choose two separate step sizes, |$\Delta _{\epsilon }$| and |$\Delta _{O}$|⁠, for emissivity and offset, respectively. For the various schemes used to tune the HMC hyper-parameters, see Hoffman & Gelman 2014, Bou-Rabee & Sanz-Serna 2017, Hoffman et al. 2021, and Sountsov & Hoffman 2021 and for their implementation, refer to the probabilistic programming frameworks such as STAN (Carpenter et al. 2017), PyMC (Salvatier et al. 2016), and pyro (Bingham et al. 2019).

A1.3 Convergence test

To quantify the convergence of the Monte Carlo sample, we adopt the Gelman–Rubin test (Gelman & Rubin 1992). To implement this test, one needs multiple Monte Carlo chains starting with sufficiently separated positions in the parameter space. One then computes the following ratio R, which should be ideally equal to one.

where V is the variance of the given parameter between the chain and W is the variance of the same parameter along the chain (Brooks & Gelman 1998; Heavens 2009). In practice, it is recommended that R should be less than 1.01 to consider the sample as converged to the distribution being sampled (Vehtari et al. 2021).

APPENDIX B: VARIABLE LEAP-FROG STEPS

To investigate the change in correlation length with Leap-Frog jumps N, we perform the HMC sampling with variable N. Following Neal (2012), we randomly draw N at each iteration from a uniform distribution, |$\mathcal {U}\left[9,15\right]$| and apply the HMC sampler on the simulated map at 353 GHz over MaskHI6 as discussed in Section 2.4. We then compute the compute the autocorrelation coefficient of emissivity |$\epsilon$| at one superpixel and the global offset O as in Section 4.2. Fig. B1 shows the autocorrelation coefficient for fixed Leap-Frog jumps |$N =10$| (Section 4.2) and variable Leap-Frog jumps |$N \in \left[9,15\right]$|⁠. It shows that for variable N, the correlation length decreases considerably for O. Variable N minimizes the correlation among the various global parameters for multiple template fit or fit with dipole contribution.