  We present an implementation of a Bayesian mixture model using Hamiltonian
Monte Carlo (HMC) techniques to search for spatial separation of Galactic dust
populations. Utilizing intensity measurements from Planck High Frequency
Instrument (HFI), we apply this model to high-latitude Galactic dust emission.
Our analysis reveals a strong preference for a spatially-varying two-population
dust model over a one-population dust model, when the latter must capture the
total variance in the sky. Each dust population is well characterized by a
single-component spectral energy distribution (SED) and accommodates small
variations. These populations could signify two distinct components, or may
originate from a one-component model with different temperatures resulting in
different SED scalings. While no spatial information is built into the
likelihood, our investigation unveils large-scale spatially coherent structures
with high significance, pointing to a physical origin for the observed spatial
variation. These results are robust to our choice of likelihood and of input
data. Furthermore, this spatially varying two-population model is the most
favored from Bayesian evidence calculations. Incorporating IRAS 100 $\mu$m to
constrain the Wein-side of the blackbody function, we find the dust populations
differ at the 2.5$\sigma$ level in the spectral index ($\beta_d$) vs.
temperature ($T_d$) plane. The presence of multiple dust populations has
implications for component separation techniques frequently employed in the
recovery of the cosmic microwave background.
