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

Debabrata Adak 1,2,3, Shabbir Shaikh4, Srijita Sinha5, Tuhin Ghosh5, Francois
Boulanger6, Guilaine Lagache7, Tarun Souradeep1,8,9 and Marc-Antoine Miville-
Deschênes10
1 Inter University Centre for Astronomy and Astrophysics, Post Bag 4,
Ganeshkhind, Pune-411007, India
2 The Institute of Mathematical Sciences, CIT Campus, Tharamani, Chennai,
Tamil Nadu 600113, India
3 Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain
4 School of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287, USA
5 National Institute of Science Education and Research, An OCC of Homi Bhabha
National Institute, Bhubaneswar 752050, Odisha, India
6 Laboratoire de Physique de l’Ecole normale supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France
7 Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France
8 Indian Institute of Science Education and Research, Dr. Homi Bhabha Road,
Ward No. 8, NCL Colony, Pashan, Pune-411008, India
9 Raman Research Institute, C. V. Raman Avenue, Bengaluru, Karnataka 560080,
India
10 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, F-91191 Gif-
sur-Yvette, France E-mail: debuadak<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a Bayesian inference method to characterise the dust emission
properties using the well-known dust-$\mathsc{Hi}$ correlation in the diffuse
interstellar medium at Planck frequencies $\nu\geq 217$ GHz. We use the
Galactic $\mathsc{Hi}$ 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 multi-
dimensional parameter space. We demonstrate that the methodology is unbiased
when applied to realistic Planck sky simulations over a 7500 deg2 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 $\mathsc{Hi}$ 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 multi-phase interstellar medium.

###### keywords:

methods: statistical, submillimetre: diffuse background, submillimetre: ISM,
cosmology: diffuse radiation

††pubyear: 2020††pagerange: Bayesian inference methodology to characterize the
dust emissivity at far-infrared and submillimeter frequencies–B

##  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 et al. 2005) are similar. While interesting in its own right,
characterisation 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 polarisation.
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
CIB is weakly polarized, at frequencies where CIB is significant, the
polarization of the CIB can be a contaminant to CMB B-mode polarization
(Lagache et al., 2020; Feng & Holder, 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 Planck111http://www.esa.int/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 obtain the 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-$\mathsc{Hi}$ correlation. Galactic dust is correlated with $\mathsc{Hi}$
21 cm line emission of neutral hydrogen at high Galactic latitude and low-
column density regions (Boulanger & Perault, 1988). Hence, the Galactic
$\mathsc{Hi}$ emission map can be used as a tracer of the Galactic dust
emission. Using Planck and IRAS data along with $\mathsc{Hi}$ 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<-25^{\circ}$).
A proper estimation of the Galactic dust emission is needed to separate the
CIB emission (see, e.g., Planck intermediate results. XLVIII. 2016; Lenz et
al. 2019; Irfan et al. 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 $\mathsc{Hi}$4PI
data (HI4PI Collaboration: et al., 2016) with the Planck intensities at 353,
545, and 857 GHz over approximately 25$\%$ 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-$\mathsc{Hi}$ correlation at high
Galactic latitude where $\mathsc{Hi}$ is the reliable Galactic dust tracer.
The underlying assumption is that the dust emission is zero where
$\mathsc{Hi}$ 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 $\mathsc{Hi}$ 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) utilises
dust-$\mathsc{Hi}$ 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-$\mathsc{Hi}$ correlation over the
patches of 15∘ diameter. Both works use the $\chi^{2}$ minimisation technique.
The efficiency of the $\chi^{2}$ minimisation 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) measures
global offset first and then estimates the Galactic dust SED parameters by
fitting the MBB spectrum to offset-corrected intensity maps. In this work, we
utilise the dust-$\mathsc{Hi}$ 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 $\mathsc{Hi}$ 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 algorithm (Duane et al., 1987). 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 et al., 2008; Jasche et al.,
2010; Jasche & Wandelt, 2013; Anderes et al., 2015), including the inference
of CMB foreground parameters (Grumitt et al., 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 summarise
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\thinspace\geq\thinspace 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 multi-frequency 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
$\mathsc{Hi}$ 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 $\mathsc{Hi}$-correlated dust emission
($I^{\mathsc{\mathsc{Hi}\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 $\texttt{H}_{2}$ gas, which is uncorrelated with the
$\mathsc{Hi}$-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 $\mathsc{Hi}$ 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 $\tt 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,

$d_{\nu}(\Omega^{j}_{i})=I^{\mathsc{\mathsc{Hi}\text{-d}}}_{\nu}(\Omega^{j}_{i})+O_{\nu}+I^{\rm
N}_{\nu}(\Omega^{j}_{i})+I^{\rm CIB}_{\nu}(\Omega^{j}_{i})+I^{\rm
R}_{\nu}(\Omega^{j}_{i}),$ (1)

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

$I^{\mathsc{\mathsc{Hi}\text{-d}}}_{\nu}(\Omega^{j}_{i})=\epsilon^{j}_{\nu}N{Hi}(\Omega^{j}_{i})\
,$ (2)

where $N{Hi}$ denotes the integrated column density of $\mathsc{Hi}$ 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 $\mathsc{Hi}$ template to trace the $\mathsc{Hi}$-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^{\mathsc{\mathsc{Hi}\text{-d}}}_{\nu}$ and $O_{\nu}$ as the signal model,
$s_{\nu}$,

$s_{\nu}(\Omega^{j}_{i})=\epsilon^{j}_{\nu}N{Hi}(\Omega^{j}_{i})+O_{\nu}\ .$
(3)

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 nonzero elements of $\boldsymbol{\Sigma_{\nu}}$ are given by

$\Sigma^{jj}_{\nu,ii^{\prime}}=\Sigma_{\nu}^{\rm
N}(\theta^{jj}_{ii^{\prime}})+\Sigma_{\nu}^{\rm
CIB}(\theta^{jj}_{ii^{\prime}})+\Sigma_{\nu}^{\rm
R}(\theta^{jj}_{ii^{\prime}}),$ (4)

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

$\mathcal{L}(\\{d_{\nu}(\Omega^{j}_{i})\\}|\\{\epsilon^{j}_{\nu},O_{\nu}\\})=\frac{1}{\sqrt{|2\pi\boldsymbol{\Sigma_{\nu}}|}}\exp{\Big{[}-\frac{\chi^{2}}{2}\Big{]}},$
(5)

where $\chi^{2}$ is

$\chi^{2}=\sum_{j}\sum_{i,i^{\prime}\subset
j}[d_{\nu}(\Omega^{j}_{i})-s_{\nu}(\Omega^{j}_{i})][\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}[d_{\nu}(\Omega^{j}_{i^{\prime}})-s_{\nu}(\Omega^{j}_{i^{\prime}})].$
(6)

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

In the model fitting, templates $N{Hi}(\Omega^{j}_{i})$ and the noise variance
$\Sigma_{\nu}$ are known, and $\\{\epsilon^{j}_{\nu},O_{\nu}\\}$ 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}(\\{\epsilon^{j}_{\nu},O_{\nu}\\}|\\{d_{\nu}(\Omega^{j}_{i})\\})$)
of the parameters given data as

$\mathcal{P}(\\{\epsilon^{j}_{\nu},O_{\nu}\\}|\\{d_{\nu}(\Omega^{j}_{i})\\})=\frac{\mathcal{L}(\\{d_{\nu}(\Omega^{j}_{i})\\}|\\{\epsilon^{j}_{\nu},O_{\nu}\\})\mathcal{P}(\\{\epsilon^{j}_{\nu},O_{\nu}\\})}{\mathcal{P}(\\{d_{\nu}(\Omega^{j}_{i})\\})},$
(7)

where $\mathcal{P}(\\{\epsilon^{j}_{\nu},O_{\nu}\\})$ is the prior probability
distribution of the parameters, and
$\mathcal{P}(\\{d_{\nu}(\Omega^{j}_{i})\\})$ is the evidence. We assume a
uniform prior for all the parameters of interest.
$\mathcal{P}(\\{d_{\nu}(\Omega^{j}_{i})\\})$ acts as a normalisation 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
$\\{\epsilon^{j}_{\nu},O_{\nu}\\}$. We have around $2\times 10^{3}$ emissivity
parameters per $N{Hi}$ 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,

$s_{\nu}(\Omega^{j}_{i})=\epsilon^{j}_{\nu}N{Hi}(\Omega^{j}_{i})+O_{\nu}+D_{\nu}(\Omega^{j}_{i}),$
(8)

where $D_{\nu}(\Omega^{j}_{i})$ 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

$D_{\nu}(\Omega^{j}_{i})=a^{\nu}_{1,0}Y_{1,0}(\Omega^{j}_{i})+a^{\nu}_{1,1}Y_{1,1}(\Omega^{j}_{i})+a^{\nu}_{1,-1}Y_{1,-1}(\Omega^{j}_{i}),$
(9)

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:

$D_{\nu}(\Omega^{j}_{i})=a^{\nu}_{1,0}Y_{1,0}(\Omega^{j}_{i})+2a^{R,\nu}_{1,1}Y^{R}_{1,1}(\Omega^{j}_{i})-2a^{I,\nu}_{1,1}Y^{I}_{1,1}(\Omega^{j}_{i}),$
(10)

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-millimeter frequency bands can also
be modelled as a linear combination of multiple $\mathsc{Hi}$ 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 $\mathsc{Hi}$ 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 $\mathsc{Hi}$ brightness temperature from the HI4PI
survey (HI4PI Collaboration: et al., 2016) to model the dust emission at
Planck HFI frequencies $\nu\geq 353$ GHz. The modelling and inference
framework used in this work can be extended to include multiple $N{Hi}$
templates with the signal modelled as

$s_{\nu}(\Omega^{j}_{i})=\sum^{N_{t}}_{t=1}\epsilon^{j,t}_{\nu}N{Hi}^{t}(\Omega^{j}_{i})+O_{\nu}+D_{\nu}(\Omega^{j}_{i}),$
(11)

where the summation is over $N_{t}$ number of $N{Hi}$ templates ($N{Hi}^{t}$),
indexed by $t$, and $\epsilon^{j,t}_{\nu}$ is corresponding emissivity.

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

###  2.3 CIB model power spectra

Figure 1: The CIB model power spectrum (in units of $\ell C_{\ell}$) obtained
by fitting the CIB model including the shot noise at four HFI frequencies
(Planck 2013 results. XXX., 2014). We use the model $C_{\ell}^{\rm CIB}$ to
compute the full covariance matrix $\Sigma_{\nu}^{\rm CIB}$.

The CIB is a relic emission from stellar-heated dust within galaxies formed
throughout cosmic history. At far-infrared/sub-millimeter 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-fit CIB model power spectra (including the shot
noise) at 217, 353, 545, and 857 GHz obtained by Planck 2013 results. XXX.
(2014). Figure 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

$\Sigma_{\nu}^{\rm
CIB}(\theta_{ii^{\prime}})=\frac{1}{4\pi}\sum(2\ell+1)\thinspace
C_{\ell,\nu\times\nu}^{\rm CIB}\thinspace B_{\ell}^{2}\thinspace
w_{\ell}^{2}\thinspace P_{\ell}(\hat{n}_{i}.\hat{n}_{i^{\prime}}),$ (12)

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 $\tt 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 $\\{\epsilon^{j}_{\nu}\\}$ and
$\\{O_{\nu}\\}$ 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 $\\{\epsilon^{j}_{\nu}\\}$
and $\\{O_{\nu}\\}$ as position variables and augments them with the momentum
variables $(p^{j,t}_{k\nu},\thinspace p_{O_{\nu}})$ to define phase space
dynamics.The Hamiltonian of the dynamics for the given problem is

$\mathcal{H}(p^{j,t}_{\nu},p_{O_{\nu}},\epsilon^{j,t}_{\nu},O_{\nu})=\frac{p^{2}_{O_{\nu}}}{2\mu_{O_{\nu}}}+\sum_{j}\frac{{p^{j,t}}^{2}_{\nu}}{2\mu^{j,t}_{\nu}}-\ln[\mathcal{P}(\\{\epsilon^{j,t}_{\nu},O_{\nu}\\})],$
(13)

where $\mathcal{P}(\\{\epsilon^{j,t}_{\nu},O_{\nu}\\})$ 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}(\\{\epsilon^{j,t}_{\nu},O_{\nu}\\})]$ is

$\ln[\mathcal{P}(\\{\epsilon^{j,t}_{\nu},O_{\nu}\\})]=-\frac{1}{2}\chi^{2}\ ,$
(14)

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. 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:

$\dot{\epsilon}^{j,t}_{\nu}=\frac{\partial\mathcal{H}}{\partial
p^{j,t}_{\nu}}=\frac{p^{j,t}_{\nu}}{\mu^{j,t}_{\nu}}\quad\text{and}\quad\dot{O}_{\nu}=\frac{\partial\mathcal{H}}{\partial
p_{O_{\nu}}}=\frac{p_{O_{\nu}}}{\mu_{O_{\nu}}},$ (15)

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:

$\dot{p}\equiv\frac{\partial\mathcal{H}}{\partial
q}=-\frac{\partial\ln(\mathcal{P})}{\partial q}.$ (16)

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

$\frac{\partial\ln(\mathcal{P})}{\partial\epsilon^{j,t}_{\nu}}=\sum_{i,i^{\prime}\subset
j}N{Hi}^{t}(\Omega^{j}_{i})[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}[d_{\nu}-s_{\nu}]^{j}_{i^{\prime}}\
,$ (17)

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

$\frac{\partial\ln(\mathcal{P})}{\partial
O_{\nu}}=\sum_{j}\sum_{i,i^{\prime}\subset
j}[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}[d_{\nu}-s_{\nu}]^{j}_{i^{\prime}}.$
(18)

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
$\\{\epsilon^{j,t}_{\nu},\thinspace O_{\nu},\thinspace a^{(R/I),\nu}_{1m}\\}$.
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

$\begin{split}\frac{\partial\ln(\mathcal{P})}{\partial
a^{(R/I),\nu}_{1m}}=(\pm)2^{m}\sum_{j}\sum_{i,i^{\prime}\subset
j}Y^{(R/I)}_{1m}(\Omega^{j}_{i})&[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}[d_{\nu}-s_{\nu}]^{j}_{i^{\prime}},\\\
&\text{for $m=0$ and $1$},\end{split}$ (19)

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. However, with the proper choice of mass matrix, $\Delta$ can be
made independent of the parameter, and one single choice of $\Delta$ suffices
(Taylor et al., 2008). This choice of $\Delta$ in the case of Gaussian
likelihoods is discussed in Taylor et al. (2008). 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 $\\{\epsilon^{j,t}_{\nu}\\}$. 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

$\mu^{j,t}_{\nu}=\sum_{i,i^{\prime}\subset
j}N{Hi}^{t}(\Omega^{j}_{i})[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}N{Hi}^{t}(\Omega^{j}_{i^{\prime}}),$
(20)

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

$\mu_{O_{\nu}}=\sum_{j}\sum_{i,i^{\prime}\subset
j}[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}},$ (21)

and the following is the mass for the dipole coefficients

$\begin{split}\mu^{(R/I)}_{a_{1m}}=2^{2m}\sum_{j}\sum_{i,i^{\prime}\subset
j}Y^{(R/I)}_{1m}(\Omega^{j}_{i})&[\Sigma_{\nu}^{-1}]^{jj}_{ii^{\prime}}Y^{(R/I)}_{1m}(\Omega^{j}_{i^{\prime}}),\\\
&\text{for $m=0$ and $1$}.\end{split}$ (22)

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.

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, $\mathsc{Hi}$ data, and the sky
mask used in the analysis. We also describe the procedure to compute the CIB
covariance matrix using the model CIB power spectrum.

###  3.1 Planck data

We use the publicly available Planck 2018 Public Release 3
(PR3222http://pla.esac.esa.int/pla) legacy intensity map at 353 GHz (Planck
2018 results. I., 2020) for our analysis. The 353 GHz intensity map has been
provided in $\tt HEALPix$333http://healpix.jpl.nasa.gov (Górski et al., 2005)
grid at $N_{\rm side}=2048$ (pixel size 1.7′) with an angular resolution of
FWHM 4.82′ (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′ (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′ (FWHM),
downgrade to $N_{\rm side}=512$ (pixel size 6.8′) resolution. We chose 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′ 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 $\mathsc{Hi}$ data

To separate the $\mathsc{Hi}$-correlated dust emission from the Planck data at
high Galactic latitude, we use the $\mathsc{Hi}$ data provided by
GASS444https://www.atnf.csiro.au/research/GASS/Data.html 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<1^{\circ}$)
within velocity range, $-400$ km s${}^{-1}<V_{\rm LSR}<500$ km s-1; where
$V_{\rm LSR}$ is the velocity of the $\mathsc{Hi}$ clouds with respect to
local standard of rest. The survey has a beam resolution of 14.5′ 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 $\mathsc{Hi}$ in the Galactic disk (or what
we call Galactic $\mathsc{Hi}$) is mixed with significant emission from the
Magellanic Stream (MS) (Nidever et al. 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 $\mathsc{Hi}$ spectra using
the velocity information (Venzmer et al., 2012). The three-dimensional model
of Kalberla & Dedes (2008) helps to distinguish the spectra associated with
the Galactic $\mathsc{Hi}$ emission and MS. Finally, the $\mathsc{Hi}$
template map is produced integrating 3D spectra over velocity range and
projected on $\tt HEALPix$ grid at $N_{\rm side}=1024$. The Galactic
$\mathsc{Hi}$ column density map $N{Hi}$ used in our analysis has an angular
resolution of 16.2′ FWHM and is projected on the $\tt HEALPix$ grid $N_{\rm
side}=512$ (pixel size 6.9′). We use the Galactic $\mathsc{Hi}$ map as a
tracer for the $\mathsc{Hi}$-correlated dust emission. The same Galactic
$\mathsc{Hi}$ 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-$\mathsc{Hi}$ correlation analysis. The total sky area of the mask is
6300 deg2 ($f_{\rm sky}=15.3\%$) where $N{Hi}<6\times 10^{20}$ cm-2, and thus
avoids the high column density regions. The unmasked sky region covers the
Galactic latitude $b\leq-25^{\circ}$. The area of 20∘ diameter centred around
($l_{MS}$, $b_{MS}$) = ($-50^{\circ}$, 0∘) is masked to avoid Magellanic
Stream (Nidever et al., 2010). Further, the bright radio sources at microwave
frequencies and infrared galaxies at 100 $\mu$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{Hi}$ 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{Hi}$ values higher than the cutoff value. We label
the mask with $N{Hi}$ 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 $\mathsc{Hi}$ 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 superpixel555Here, the threshold is 85 unmasked
subpixels within superpixels at $N_{\rm side}=32$.. 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:

  1. 1.

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) $\tt HEALPix$ grid. This map is obtained through the
dust-$\mathsc{Hi}$ correlation analysis over 15∘ circular patches in diameter
centred on $\tt 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.

  2. 2.

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

$\epsilon_{\nu}(\Omega^{j}_{i})=\epsilon_{353}(\Omega^{j}_{i})\left(\frac{\nu}{353}\right)^{\beta_{\rm
d}}\frac{B_{\nu}({T_{\rm d}})}{B_{353}({T_{\rm d}})},$ (23)

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

  3. 3.

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 $\mathsc{Hi}$ template.
We add the global offset values from Planck intermediate results. XLVIII.
(2016) to produce the signal model maps at all HFI frequencies.

  4. 4.

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 CIB is correlated between two frequencies, analysing
individual frequency maps entails neglecting the correlation between the CIB
emission at different frequencies.

  5. 5.

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

  6. 6.

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.

Figure 2: The auto-correlation coefficient, $\rho(t)$, at 353 GHz for the samples of the emissivity in one superpixel ($\epsilon$) and the global offset ($O$) as a function of the lag $(t)$ for two different choices of Leap-Frog step size ($\Delta_{O}$) for the global offset and fixed $\Delta_{\epsilon}=0.1$. For $\Delta_{O}=0.1$, the correlation length for $\epsilon$ and $O$ are $8$ and $220$, respectively. For $\Delta_{O}=1.0$, the correlation length for $\epsilon$ in the same pixel and $O$ are $10$ and $102$, respectively. Table 1: Results of analysis of simulated maps at four Planck HFI frequencies. The global offset values with corresponding 1$\sigma$ error bars are estimated over five sky masks, tracing low to intermediate $\mathsc{Hi}$ column density regions. The inference of offset is stable with respect to sky coverage, with a slight decrease in the uncertainty with higher $f_{\rm sky}$.  Frequency | Input offsets | Recovered Offsets [ kJy sr-1]
---|---|---
$[$ GHz$]$ | [ kJy sr-1] |
|  | Sky masks
|  | MaskHI2 | MaskHI3 | MaskHI4 | MaskHI5 | MaskHI6
|  |  $f_{\rm sky}$ [%]
|  | 7.3 | 11.5 | 13.9 | 15.0 | 15.3
217 | 40 |  $40.0\pm 0.1$ |  $40.0\pm 0.1$ |  $40.0\pm 0.1$ |  $40.0\pm 0.1$ |  $40.0\pm 0.1$
353 | 120 | $119.6\pm 0.4$ | $120.1\pm 0.4$ | $119.9\pm 0.3$ | $119.9\pm 0.3$ | $120.0\pm 0.2$
545 | 330 | $330.9\pm 0.9$ | $330.5\pm 0.7$ | $330.3\pm 0.7$ | $330.6\pm 0.7$ | $330.6\pm 0.6$
857 | 550 | $550.5\pm 1.2$ | $550.6\pm 1.0$ | $550.6\pm 0.9$ | $550.7\pm 0.9$ | $550.6\pm 0.9$

###  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 summarise 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.

Figure 3: The joint and marginalised probability distributions of dust
emissivities (expressed in kJy sr-1$(10^{20}\thinspace\text{cm}^{-2})^{-1}$)
at three representative superpixels (pixel index as super-script) and the
global offset (in kJy sr-1) at 353 GHz. The red lines mark the posterior mean,
and the orange lines depict the input values of the respective parameters. The
contours mark the $68\%$ and $90\%$ regions of the joint distributions. The
vertical black dashed lines in the histogram mark $16$, $50$, and $84$
percentiles of the distribution.

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. Figure 2
presents the auto-correlation 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
$\\{\Delta_{\epsilon},\Delta_{O}\\}$. When step size for $\epsilon$ and ${O}$
are the same, $\\{\Delta_{\epsilon},\Delta_{O}\\}=\\{0.1,0.1\\}$, 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 $\\{\Delta_{\epsilon},\Delta_{O}\\}=\\{0.1,1.0\\}$, 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
$\\{\Delta_{\epsilon},\Delta_{O}\\}$, 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 A.1). 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.

Figure 4: The histograms of the normalised difference of input emissivities
and the posterior mean emissivities,
$\delta^{j}=(\epsilon^{j}_{inp}-\epsilon^{j}_{out})/\sigma_{\epsilon^{j}_{out}}$
for four Planck frequencies from 217 to 857 GHz. Vertical lines depict the
normalised difference between the input and output global offset at respective
Planck frequencies.  Figure 5: The recovered dust emissivity
$\epsilon^{j}_{\nu}$ and its standard deviation $\sigma_{\epsilon^{j}_{\nu}}$
per superpixel ($j$) as a function of mean Galactic $\mathsc{Hi}$ column
density ($\langle N{Hi}(\Omega^{j}_{i})\rangle$). The mean $\langle..\rangle$
is taken over all unmasked subpixels ($i$) that fall within a given
superpixel. The “dot" symbol represents $\epsilon^{j}_{\nu}$ and “plus"
represents $\sigma_{\epsilon^{j}_{\nu}}$ at all Planck frequencies considered
in this analysis. Figure 6: The left panel shows the input simulated map at
353 GHz (in units of kJy sr-1) over MaskHI6. The middle panel shows the signal
model map (in kJy sr-1) derived from the mean values of the emissivity and the
offset obtained using the HMC sampler. The right panel depicts the residual
map (in kJy sr-1) obtained by subtracting the signal model map from the input
map at 353 GHz. The residual map has contributions from the CIB and the
instrumental noise.  Figure 7: Figure shows the histograms of the residual map
(black curve) after subtracting the signal model from the input simulated map
at Planck HFI frequencies over the unmasked pixels of eMaskHI6. The expected
histograms of the residuals (magenta curve) are produced from a single
realization of instrument noise and the CIB at the respective frequency over
the same sky mask. The agreement between the observed and expected residual
distribution validates that our HMC methodology is unbiased.

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 Figure 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 marginalised
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.

Figure 8: The ratio of $1\sigma$ uncertainty in dust emissivity with and
without the CIB contribution in the noise covariance term. The ratio is
plotted for all the superpixels over the unmasked sky area against their mean
$N{Hi}$ value (“plus" marker). The horizontal lines represent the ratio of the
uncertainty in the offset parameter, which is a single number at a given
frequency, with and without the CIB in the noise covariance.

To quantify the goodness of fit, we compute the reduced $\chi^{2}$ considering
the pixels included in the analysis, which are a subset of those in the
MaskHIQ mask. This is due to excluding superpixels that do not satisfy the
threshold criteria on 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 find the
$\chi^{2}$ per degree of freedom at best-fit parameters equal to 1.00, 1.02,
1.06, and 1.19 at respective Planck HFI frequencies over the eMaskHI6 mask.

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 Figure 4, we present the distribution of the standard deviation
($\sigma_{\epsilon^{j}_{out}}$) normalised 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-fit values of output dust emissivities for all superpixels are unbiased.

In Figure 5, we show the mean and standard deviation of dust emissivity per
superpixel obtained from the dust-$\mathsc{Hi}$ correlation analysis as a
function of mean Galactic $\mathsc{Hi}$ column density at those superpixels.
The uncertainty decreases with an increase in mean Galactic $\mathsc{Hi}$
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, which is equal to $1/\mu_{O_{\nu}}$ (see Equation 21), is
inversely proportional to $N{Hi}$.

Owing to the linear nature of the signal model, the signal that corresponds to
the best-fit values of emissivity and the offset is also the best-fit signal
corresponding to the maximum of the posterior. We obtain the best-fit
intensity map corresponding to the signal model using the mean values of the
dust emissivity and the global offset following Equation 3. In Figure 6, we
show the simulated map at 353 GHz along with maps of the best-fit model and
the residual. The residual map is the difference between the input and best-
fit 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
Figure 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 CIB noise in our analysis, we redo the analysis
without the CIB contribution in the noise covariance matrix. In Figure 8, we
compare the standard deviation of the dust emissivity and offset inferred with
and without considering CIB in the noise covariance matrix. We plot the ratio
of these two standard deviations for all the superpixels against the mean
$\langle N{Hi}\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 CIB as a
function of $N{Hi}$ 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 Figure 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 $\mathsc{Hi}$-correlated dust
emission and Galactic residuals, one can expect it to dominate over 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 $\mathsc{Hi}$ 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{Hi}$ templates. We use
the MS and Galactic $\mathsc{Hi}$ 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^{\mathsc{Hi}}_{353}\rangle$,
where $\langle\epsilon^{\mathsc{Hi}}_{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^{\mathsc{Hi}}_{353}\rangle=37.3\thinspace\thinspace$kJy
sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ and
$\epsilon^{\mathrm{MS}}_{353}=0.37\thinspace$kJy
sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$, respectively, and the offset is
$120\thinspace$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 $\mathsc{Hi}$ template and the MS template over MaskHI6 performing
the HMC methodology as discussed in Section 2.4. Figure 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\thinspace\thinspace$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\thinspace\thinspace$kJy
sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ and $7.2\thinspace\thinspace$kJy
sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ for the Galactic $\mathsc{Hi}$ template
while $0.39\thinspace\thinspace$kJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ and
$0.27\thinspace\thinspace$kJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ for the MS
template.

Figure 9: The histograms of the normalised difference of input and the
posterior mean emissivities for Galactic $\mathsc{Hi}$ and MS template at 353
GHz. The vertical line shows the normalised difference between the input and
output global offset at the same frequency.

###  4.4 Validation with a residual CMB dipole term

Table 2: Recovered global offset ($O$), dipole amplitude ($A_{\mathrm{dip}}$) and Galactic longitude ($l_{\mathrm{dip}}$) and latitude ($b_{\mathrm{dip}}$) estimated over five different $N{Hi}$ cutoff sky masks for simulated map at 217 GHz. The corresponding input values are respectively $O=40\thinspace\thinspace$kJy sr-1$,A_{\mathrm{dip}}=3.9\thinspace\thinspace$kJy sr-1$,l_{\mathrm{dip}}=264.0\degree,b_{\mathrm{dip}}=48.3\degree$. Column $\delta$ represents the difference between the input and the recovered parameter value in units of $1\sigma$ uncertainty. Sky masks |  $O\thinspace\left[\mathrm{kJy\thinspace sr^{-1}}\right]$ | $A_{\mathrm{dip}}\thinspace\left[\mathrm{kJy\thinspace sr^{-1}}\right]$ | $l_{\mathrm{dip}}\thinspace\left[\mathrm{degree}\right]$ | $b_{\mathrm{dip}}\left[\mathrm{degree}\right]$
---|---|---|---|---
|  Recovered | $\delta$ | Recovered | $\delta$ | Recovered | $\delta$ | Recovered | $\delta$
MaskHI2 | $39.7\pm 1.6$ | $0.2$ | $4.1\pm 0.8$ | $-0.3$ | $287.2\pm 9.4$ | $-2.5$ | $38.4\pm 24.6$ | $0.4$
MaskHI3 | $39.3\pm 0.9$ | $0.7$ | $3.6\pm 0.5$ | $0.6$ | $271.2\pm 6.5$ | $-1.1$ | $34.0\pm 15.9$ | $0.9$
MaskHI4 | $39.2\pm 0.7$ | $1.1$ | $3.8\pm 0.4$ | $0.2$ | $273.0\pm 4.2$ | $-2.1$ | $30.5\pm 12.5$ | $1.4$
MaskHI5 | $39.8\pm 0.6$ | $0.4$ | $3.9\pm 0.4$ | $-0.1$ | $273.1\pm 4.7$ | $-1.9$ | $39.7\pm 10.0$ | $0.9$
MaskHI6 | $39.2\pm 0.6$ | $1.2$ | $3.7\pm 0.3$ | $0.6$ | $271.4\pm 4.2$ | $-1.7$ | $31.0\pm 10.9$ | $1.6$

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\thinspace$kJy sr-1,
corresponding to $8\thinspace\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\degree,48.3\degree\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\thinspace\thinspace$kJy sr-1 (same as in Table 1), and
the input dust emissivity averaged over the superpixels for 217 GHz is
$8.15\thinspace$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\thinspace$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 normalised histogram,
indicating the method is unbiased. The harmonic coefficients are related to
the dipole amplitude and Galactic longitude and latitude, respectively, as

$\displaystyle A_{\mathrm{dip}}$
$\displaystyle=\sqrt{\frac{3}{4\pi}}\sqrt{a_{1,0}+2a^{R}_{1,1}+2a^{I}_{1,1}},$
$\displaystyle l_{\mathrm{dip}}$
$\displaystyle=\arctan\left(-\frac{a^{I}_{1,1}}{a^{R}_{1,1}}\right),$ (24)
$\displaystyle b_{\mathrm{dip}}$
$\displaystyle=\frac{\pi}{2}-\arccos\left(\frac{a_{1,0}}{A_{\mathrm{dip}}}\sqrt{\frac{3}{4\pi}}\right).$

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 marginalised
probability distributions of the global offset, dipole amplitude and direction
are shown in Figure 10. The corresponding mean values with standard deviations
are obtained as $A_{\mathrm{dip}}=3.7\pm 0.3\thinspace\thinspace$kJy sr-1 and
direction $\left(l_{\mathrm{dip}},b_{\mathrm{dip}}\right)=\left(271.4\pm
4.2\degree,31.0\pm 10.9\degree\right)$.

Figure 10: The joint and marginalised probability distributions of global
offset and the three dipole parameters at 217 GHz. The global offset and
residual CMB dipole amplitude are expressed in kJy sr-1 unit and the residual
CMB dipole direction in degrees. Details are the same as Figure 3.

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 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 2 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\thinspace\thinspace$kJy sr-1 (in intensity units). We
performed the same analysis and recovered the offset $O=119.0\pm
2.4\thinspace\thinspace$kJy sr-1 and dipole amplitude $A_{\mathrm{dip}}=3.7\pm
1.4\thinspace\thinspace$kJy sr-1 and direction as
$\left(l_{\mathrm{dip}},b_{\mathrm{dip}}\right)=\left(277.4\pm
31.1\degree,14.9\pm 43.2\degree\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 Figure 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. We do not consider the
residual CMB dipole contribution in the signal model for this analysis. 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. The choices
made for the HMC parameters, such as Leap-Frog step sizes and the number of
jumps in each Leap-Frog loop while analysing simulated maps, are also suitable
for analyzing the actual 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.

Figure 11: The mean (_left_) and the standard deviation (_right_) of the dust
emissivity obtained from the Planck 353 GHz intensity map. Both the maps are
in MJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ units. These are for the analysis
with the MaskHI6 mask, which has $f_{\rm sky}=15.3\%$. Figure 12: _Upper
panel_ : Histogram of mean dust emissivity obtained in the Planck 353 GHz
intensity map analysis over MaskHI6. The map of these dust emissivity values
is given in the _left panel_ of Figure 11. _Lower panel_ : Variation of dust
emissivity as a function of mean dust column density. $\langle N{Hi}\rangle$
represents the bin average $\mathsc{Hi}$ column density and
$\langle\epsilon\rangle$ are the average over the superpixels within the same
$N{Hi}$ bins. The error bars are the standard deviations in the respective
bins.

In Figure 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{Hi}$ 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 Figure 11 is $0.031\thinspace\thinspace$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\thinspace\thinspace$MJy
sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$. Note that, as shown in Figure 12, the
spatial distribution of the dust emissivity is slightly skewed towards higher
values. The upper panel of Figure 12 shows the histogram of dust emissivity
values in the $\epsilon$ map depicted in Figure 11. The lower panel of Figure
12 shows the variation of dust emissivity as a function of $N{Hi}$. 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{Hi}$ bin.
The average dust emissivity increases with increasing $\mathsc{Hi}$ column
density. Our spatial average value of the dust emissivity is $\approx 20\%$
lower than the value quoted in Planck intermediate results. XVII. (2014). This
difference could arise due to (1) different aperture size considered for the
dust-$\mathsc{Hi}$ 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\thinspace\thinspace$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 $\tt 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 Figure 13. The leftmost panel in Figure 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 Figure
13 shows the difference between the CMB-subtracted Planck 353 GHz intensity
map and signal model map over the same sky mask.

Figure 13: Results of analysis with Planck data at 353 GHz in the MJy sr-1
unit. _First panel_ : the CMB subtracted Planck intensity map at 353 GHz over
the South Galactic Pole region, _Second panel_ : the estimate of the signal
map obtained using the posterior mean emissivity and the offset, _Third panel_
: the residual map, which is the difference between the difference map and the
signal estimate, _Fourth panel_ : histogram of the residual map for the
analysis with different sky masks. The dashed blue curve shows the expected
residual based on the CIB and the instrument noise.

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-fit 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 Figure 7). The fourth panel of Figure 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 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 localised pixels in
which dust is uncorrelated with the Galactic $\mathsc{Hi}$.

Table 3 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 $\mathsc{Hi}$ 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{Hi}$ 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 $\mathsc{Hi}$ column density map,
(2) the Planck Galactic zero-level is estimated over a larger sky fraction
(28% of the sky) than considered here, and (3) the difference in the local
velocity cloud $\mathsc{Hi}$ column density threshold ($<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
$\mathsc{Hi}$ and MS templates using the methodology discussed in Section
2.2.2. We infer the offset $O=0.1281\pm 0.0003\thinspace$MJy sr-1 at 353 GHz
over MaskHI6 for the two-template fit. The average of best-fit values of
recovered emissivities for Galactic $\mathsc{Hi}$ and MS templates are
$0.031\thinspace\thinspace$MJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ and
$-0.6\times 10^{-4}\thinspace$MJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$,
respectively. The corresponding $1\sigma$ deviation on them are
$0.007\thinspace\thinspace$MJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$ and
$1.8\times 10^{-4}\thinspace$MJy sr-1$(10^{20}\mathrm{cm^{-2}})^{-1}$,
respectively.

Table 3: Results of analysis of Planck intensity map at 353 GHz. We list the global offset ($O$) values and the mean dust emissivity over the superpixels ($\langle\epsilon^{j}\rangle$) estimated over different sky masks corresponding to different $N{Hi}$ cutoffs. Column $\sigma$ represents the uncertainty on the mean dust emissivity $\langle\epsilon^{j}\rangle$. We also include the estimates from the Planck collaboration analysis for comparison. Sky masks ($f_{\rm{sky}}\left[\%\right]$) |  $O$ $\left[\mathrm{MJy\thinspace sr^{-1}}\right]$ |  $\langle\epsilon^{j}\rangle$ $\left[\mathrm{MJy\thinspace sr^{-1}}(10^{20}\mathrm{cm^{-2}})^{-1}\right]$ |  $\sigma$ $\left[\mathrm{MJy\thinspace sr^{-1}}(10^{20}\mathrm{cm^{-2}})^{-1}\right]$
---|---|---|---
MaskHI2 (7.3) | $0.1358\pm 0.0005$ | $0.023$ | $0.00022$
MaskHI3 (11.5) | $0.1332\pm 0.0004$ | $0.026$ | $0.00019$
MaskHI4 (13.9) | $0.1311\pm 0.0004$ | $0.029$ | $0.00018$
MaskHI5 (15.0) | $0.1294\pm 0.0003$ | $0.030$ | $0.00017$
MaskHI6 (15.3) | $0.1284\pm 0.0003$ | $0.031$ | $0.00016$
Planck analysis (Planck intermediate results. XVII., 2014) (Planck 2018 results. I., 2020) | $0.130\pm 0.038$ | $0.039$ | $0.0014$

##  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 utilise the
dust-$\mathsc{Hi}$ correlation to model the Galactic thermal dust emission,
using the integrated Galactic $\mathsc{Hi}$ 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 7500
deg2 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 is 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, the uncertainty that 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
$\mathsc{Hi}$ 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 $\mathsc{Hi}$ 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
Figure 11, Figure 12, Figure 13, and Table 3. As shown in Figure 11, we infer
the emissivity of the dust over the sky area $\approx 7500$ deg2 and its
variation at scales larger than 1.8∘. For the region of interest over $7500$
deg2 area centered around the southern Galactic pole with $N{Hi}$ threshold
$N{Hi}<6\times 10^{20}\text{cm}^{-2}$, the spatial average of dust emissivity
is $0.031\leavevmode\nobreak\ \mathrm{MJy\thinspace
sr^{-1}}(10^{20}\mathrm{cm^{-2}})^{-1}$ with standard deviation of
$0.007\leavevmode\nobreak\ \mathrm{MJy\thinspace
sr^{-1}}(10^{20}\mathrm{cm^{-2}})^{-1}$. The inferred offset for MaskHI6 is
$0.1284\pm 0.0003\thinspace\thinspace$MJy sr-1 is close to the monopole of the
Béthermin et al. (2012) CIB model value $0.130\pm
0.038\thinspace\thinspace$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{Hi}$ 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{Hi}$) does depend on $N{Hi}$ 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-$\mathsc{Hi}$ 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 $\mathsc{Hi}$ 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 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.

## Acknowledgements

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 $\mathsc{Hi}$ 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 $\tt 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. D. Adak acknowledges UGC for
providing Ph.D. fellowship when a major part of this project is done. D. Adak
acknowledges IMSc for providing postdoc fellowship for one year when the rest
of the project is completed. S. Shaikh 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. S.
Sinha would like to thank NISER Bhubaneswar for the postdoctoral fellowship.
G. Lagache 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 $\mathsc{Hi}$
column density map is downloaded from https://www.astro.uni-
bonn.de/hisurvey/index.php. The different sky masks used in the analysis are
publicly available at https://github.com/shabbir137/dust_emissivity.git. We
will make the dust model map and the residual map at 353 GHz available after
publication at the same location.

## References

  * Adak et al. (2020) Adak D., Ghosh T., Boulanger F., Haud U., Kalberla P., Martin P. G., Bracco A., Souradeep T., 2020, A&A, 640, A100
  * Ade et al. (2014) Ade P. A. R., et al., 2014, Astron. Astrophys., 571, A18
  * Ade et al. (2016) Ade P. A. R., et al., 2016, Astron. Astrophys., 594, A23
  * Anderes et al. (2015) Anderes E., Wandelt B., Lavaux G., 2015, Astrophys. J., 808, 152
  * Béthermin et al. (2012) Béthermin M., et al., 2012, The Astrophysical Journal, 757, L23
  * Boulanger & Perault (1988) Boulanger F., Perault M., 1988, ApJ, 330, 964
  * D’Onghia & Fox (2016) D’Onghia E., Fox A. J., 2016, Annual Review of Astronomy and Astrophysics, 54, 363
  * Duane et al. (1987) Duane S., Kennedy A. D., Pendleton B. J., Roweth D., 1987, Phys. Lett., B195, 216
  * Feng & Holder (2020) Feng C., Holder G., 2020, Astrophys. J., 897, 140
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statist. Sci., 7, 457
  * Ghosh et al. (2017) Ghosh T., et al., 2017, A&A, 601, A71
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wand elt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Grumitt et al. (2020) Grumitt R. D. P., Jew L. R. P., Dickinson C., 2020, Mon. Not. Roy. Astron. Soc., 496, 4383
  * HI4PI Collaboration: et al. (2016) HI4PI Collaboration: et al., 2016, A&A, 594, A116
  * Hajian (2007) Hajian A., 2007, Phys. Rev. D, 75, 083525
  * Heavens (2009) Heavens A., 2009, arXiv e-prints, p. arXiv:0906.0664
  * Irfan et al. (2019) Irfan M. O., Bobin J., Miville-Deschênes M.-A., Grenier I., 2019, A&A, 623, A21
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, Astrophys. J., 779, 15
  * Jasche et al. (2010) Jasche J., Kitaura F. S., Wandelt B. D., Ensslin T. A., 2010, Mon. Not. Roy. Astron. Soc., 406, 60
  * Jego et al. (2023) Jego B., Alonso D., García-García C., Ruiz-Zapatero J., 2023, MNRAS, 520, 583
  * Kalberla & Dedes (2008) Kalberla P. M. W., Dedes L., 2008, A&A, 487, 951
  * Kalberla et al. (2005) Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M., Bajaja E., Morras R., Pöppel W. G. L., 2005, A&A, 440, 775
  * Kalberla et al. (2010) Kalberla P. M. W., et al., 2010, A&A, 521, A17
  * Lagache et al. (2005) Lagache G., Puget J.-L., Dole H., 2005, Ann. Rev. Astron. Astrophys., 43, 727
  * Lagache et al. (2020) Lagache G., Béthermin M., Montier L., Serra P., Tucci M., 2020, Astron. Astrophys., 642, A232
  * Larsen et al. (2016) Larsen P., Challinor A., Sherwin B. D., Mak D., 2016, Phys. Rev. Lett., 117, 151102
  * Lenz et al. (2019) Lenz D., Doré O., Lagache G., 2019, ApJ, 883, 75
  * Maniyar et al. (2019) Maniyar A. S., Lagache G., Béthermin M., Ilić S., 2019, Astron. Astrophys., 621, A32
  * McClure-Griffiths et al. (2009) McClure-Griffiths N. M., et al., 2009, ApJS, 181, 398
  * Neal (2012) Neal R. M., 2012, arXiv e-prints, p. arXiv:1206.1901
  * Nidever et al. (2008) Nidever D. L., Majewski S. R., Butler Burton W., 2008, ApJ, 679, 432
  * Nidever et al. (2010) Nidever D. L., Majewski S. R., Butler Burton W., Nigra L., 2010, ApJ, 723, 1618
  * Planck 2013 results. IX. (2014) Planck 2013 results. IX. 2014, A&A, 571, A9
  * Planck 2013 results. VIII. (2014) Planck 2013 results. VIII. 2014, A&A, 571, A8
  * Planck 2013 results. XI. (2014) Planck 2013 results. XI. 2014, A&A, 571, A11
  * Planck 2013 results. XXVII. (2014) Planck 2013 results. XXVII. 2014, A&A, 571, A27
  * Planck 2013 results. XXX. (2014) Planck 2013 results. XXX. 2014, A&A, 571, A30
  * Planck 2015 results. VIII. (2016) Planck 2015 results. VIII. 2016, A&A, 594, A8
  * Planck 2015 results. X (2016) Planck 2015 results. X 2016, A&A, 594, A10
  * Planck 2018 results. I. (2020) Planck 2018 results. I. 2020, A&A, 641, A1
  * Planck 2018 results. III. (2020) Planck 2018 results. III. 2020, A&A, 641, A3
  * Planck 2018 results. IV. (2020) Planck 2018 results. IV. 2020, A&A, 641, A4
  * Planck 2018 results. XI. (2020) Planck 2018 results. XI. 2020, A&A, 641, A11
  * Planck early results. XVIII. (2011) Planck early results. XVIII. 2011, A&A, 536, A18
  * Planck early results. XXIV. (2011) Planck early results. XXIV. 2011, A&A, 536, A24
  * Planck intermediate results. LVII. (2020) Planck intermediate results. LVII. 2020, A&A, 643, A42
  * Planck intermediate results. XLVIII. (2016) Planck intermediate results. XLVIII. 2016, A&A, 596, A109
  * Planck intermediate results. XVII. (2014) Planck intermediate results. XVII. 2014, A&A, 566, A55
  * Planck intermediate results. XXII (2015) Planck intermediate results. XXII 2015, A&A, 576, A107
  * Puget et al. (1996) Puget J. L., Abergel A., Bernard J. P., Boulanger F., Burton W. B., Desert F. X., Hartmann D., 1996, A&A, 308, L5
  * Taylor et al. (2008) Taylor J. F., Ashdown M. A. J., Hobson M. P., 2008, Mon. Not. Roy. Astron. Soc., 389, 1284
  * Venzmer et al. (2012) Venzmer M. S., Kerp J., Kalberla P. M. W., 2012, A&A, 547, A12
  * van Engelen et al. (2015) van Engelen A., et al., 2015, Astrophys. J., 808, 7

##  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 the Hamiltonian dynamics to act as the
process of generating the proposed point. 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 parameters obey the Gaussian probability distribution with covariance
determined 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}(\\{q_{i}\\})$, the probability distribution that is
sampled in HMC is the joint probability distribution of ${p_{i}}$ and
${q_{i}}$:

$\displaystyle\mathcal{P}(\\{q_{i},p_{i}\\})$ $\displaystyle\equiv$
$\displaystyle\frac{\exp(-\mathcal{H})}{\sqrt{|2\pi\boldsymbol{\mu}|}}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{|2\pi\boldsymbol{\mu}|}}\exp\Big{[}-\frac{\textbf{p}^{T}\boldsymbol{\mu}^{-1}\textbf{p}}{2}\Big{]}\mathcal{P}(\\{q_{i}\\}).$

Here, p is the vector of momenta $\\{p_{i}\\}$ associated with the parameters
and $\mathcal{H}$ is the Hamiltonian comprising of a Kinetic Energy term and a
Potential Energy term:

$\mathcal{H}(\\{q_{i},p_{i}\\})=\frac{1}{2}\textbf{p}^{T}\boldsymbol{\mu}^{-1}\textbf{p}-\ln[\mathcal{P}(\\{q_{i}\\})]\
.$ (26)

In the above equation, the first term on the right-hand side is the kinetic
energy term, and the second term, $-\ln[\mathcal{P}(\\{q_{i}\\})]$, acts as
potential energy. $\boldsymbol{\mu}$ is the mass matrix corresponding to the
set of parameter $\\{q_{i}\\}$. 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

$\mathcal{H}(\\{q_{i},p_{i}\\})=\sum_{i}\frac{p_{i}^{2}}{2\mu_{i}}-\ln[\mathcal{P}(\\{q_{i}\\})]\
.$ (27)

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

$\dot{q}\equiv\frac{dq}{dt}=\frac{\partial\mathcal{H}}{\partial
p}\quad\text{and}\quad\dot{p}\equiv\frac{dp}{dt}=-\frac{\partial\mathcal{H}}{\partial
q}.$ (28)

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

$p(t+\Delta/2)=p(t)+\frac{\Delta}{2}\dot{p}(q(t));$ (29)
$q(t+\Delta)=q(t)+p(t+\Delta/2)\Delta/\mu;$ (30)
$p(t+\Delta)=p(t+\Delta/2)+\frac{\Delta}{2}\dot{p}(q(t+\Delta)).$ (31)

In the above equations, $\dot{p}(q(t))$ is the momentum derivative substituted
from Equation 17 or Equation 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 $\\{\epsilon^{j,t}_{\nu},\thinspace O_{\nu},\thinspace
a^{(R/I),\nu}_{1m}\\}$) and $p_{i}$.

1:Initialise $q^{(0)}$

2:for $k=1...N_{s}$ do

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

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

5: for $j=1...N$ do

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

7: end for

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

9: Accept $q^{*},p^{*}$ with the probability
$\rm{min}(1,e^{-H(q^{*},p^{*})+H(q^{(k-1)},p^{(k-1)})})$. 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 $\\{q^{(0)}_{i}\\}$. The corresponding momenta
$\\{p^{k}_{i}\\}$ are drawn from a multi-normal 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. Hence, each proposed point in phase space need not be exactly on
the same phase space trajectory on which the initial point lay. To mitigate
this numerical artefact, 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 marginalisation 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})$.

###  A.1 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.

##### Burn-in:

We neglect some of the initial samples from the chain as Burn-in. To
quantitatively determine the Burn-in sample’s size, we compute the $\chi^{2}$
per degree of freedom ($d.o.f$) using a model map estimated with the
increasing number of samples in the chains of the parameters. In each bin, the
model map has been built from the posterior mean of the parameter chain with
the increasing number of samples.

##### 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,

$L=1+2\sum_{t=1}^{t_{\rm max}}\rho(t),$ (32)

where $\rho(t)$ is the auto-correlation 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.

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.

##### 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.

$R=\frac{V}{W},$ (33)

where $V$ is the variance of the given parameter between the chain and $W$ is
the variance of the same parameter along the chain. In practice, it is
recommended that $R$ should be less than 1.2 to consider the sample as
converged to the distribution being sampled (Heavens, 2009).

##  Appendix B Variable Leap-Frog steps

Figure 14: Same as Figure 2, but for fixed and variable values of Leap-Frog
jumps $N$. For a fixed $N=10$, the correlation length for $\epsilon$ and $O$
are $10$ and $102$, respectively. For variable $N\in\left[9,15\right]$, the
correlation length for $\epsilon$ and $O$ are $7$ and $57$, respectively.

To investigate the change in correlation length with Leap-Frog jumps $N$, we
perform the HMC sampling with variable $N$. 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 auto-correlation
coefficient of emissivity $\epsilon$ at one superpixel and the global offset
$O$ as in Section 4.2. Figure 14 shows the auto-correlation 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$ minimises the correlation
among the various global parameters for multiple template fit or fit with
dipole contribution.