title: AzTEC Survey of the Central Molecular Zone: Increasing Spectral Index of Dust with Density
authors: Yuping Tang Q. Daniel Wang Grant W. Wilson
text: 
  
  
  
  
  
  
  
  Department of Astronomy, University of Massachusetts Amherst, Amherst, 01002, USA
  
  
  
  
  
  
  
  
  The Central Molecular Zone (CMZ) of our Galaxy hosts an extreme environment analogous to that found in typical starburst galaxies in the distant universe. In order to understand dust properties in environments like our CMZ, we present results from a joint SED analysis of our AzTEC/Large Millimeter Telescope survey, together with existing Herschel far-IR data on the CMZ, from a wavelength range of 160 μ m to 1.1 mm. We include global foreground and background contributions in a novel Bayesian modeling that incorporates the Point Spread Functions (PSFs) of the different maps, which enables the full utilization of our high resolution (10.5”) map at 1.1 mm and reveals unprecedentedly detailed information on the spatial distribution of dusty gas across the CMZ. There is a remarkable trend of increasing dust spectral index β, from 2.0-2.4, toward dense peaks in the CMZ, indicating a deficiency of large grains or a fundamental change in dust optical properties. This environmental dependence of β could have a significant impact on the determination of dust temperature in other studies. Depending on how the optical properties of dust deviate from the conventional model, dust temperatures could be underestimated by 10-50% in particularly dense regions. 
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
   
  
  Observations have revealed that the central ∼ 200 parsec region, or the main body of the so-called Central Molecular Zone (CMZ) of our Galaxy has an extreme gaseous environment, which may be common to the nuclear regions of many galaxies <cit.>. The CMZ is characterized by dense (n_H_2≳ 10^4 cm^-3), warm (T ≈ 60-100 K) <cit.> molecular gas with violent turbulent motions <cit.>. The magnetic fields <cit.> and the flux density of cosmic rays <cit.> in the CMZ are larger here than anywhere else in the the Galactic disk. As a result, the CMZ hosts an environment with conditions similar to those observed in high redshift starburst galaxies <cit.>. 
  
  
  The tight connections between gas conditions in the CMZ and high-redshift starburst galaxies highlight it as a template for verification/calibration of dust models in extreme environments. In modeling of high-redshift starburst galaxies <cit.>, optical properties of dust grains are conventionally adopted from those inferred in the local environments. The dust absorption curve is normally simplified as a single power-law from far infrared (FIR) to submillimeter wavelengths, characterized by a spectral index β=1.5-2. Our ability to constrain dust properties in distant starburst galaxies is limited by the lack of spatial resolution <cit.>. As a matter of fact, even in the local universe, studies of the spectral energy distributions (SEDs) of dust emission 
  have not clearly established how the optical properties of dust vary in different environments. On small scales, β≈ 1 is observed in proto-planetary and proto-stellar disks <cit.> and is commonly attributed to ≳ 1 mm size large grains. In dense molecular clouds and the diffuse ISM, a wide variety of β is observed, from 0.8 to >2 <cit.>. The origin of this diversity is debated.
  So far, observations suggest an anti-correlation between dust temperature and β or a positive correlation between gas density n_H_2 and β at long wavelengths λ≳ 500 μ m,  over the range from the diffuse ISM to cold dense clumps <cit.>. At short wavelengths ≲ 200-500 μ m, however, an inverse trend is observed <cit.>,  i.e., a flattening of the dust absorption curve toward dense regions. 
  It has been suggested that radiative transfer effects <cit.> and parameter degeneracies <cit.> could be responsible. The wavelength dependent change of β is intriguing, and cannot be reproduced by classic models of dust growth <cit.>, which predict a negative n_H2-β correlation extending to millimeter wavelengths. Recently, two new models have been proposed to solve this problem: a) accretion of small hydrogenated amorphous carbon onto large grains with updated optical properties of the hydrogenated amorphous carbons <cit.> and b) an intrinsic dependency of the dust absorption curve on the dust temperature <cit.>. These two models could be potentially distinguishable from observations, as the first scenario suggests a density dependency of β, and the second scenario suggests a temperature dependency. Nevertheless, observations with wide coverage in the n_H_2-T_dust plane are required. It is also possible that turbulence <cit.> is a factor affecting dust properties, which could potentially enhance shattering of large dust grains or suppress coagulation of small grains, especially in systems like the CMZ and high-redshift starburst galaxies where gas motions are extreme.  
  
  Submillimeter/millimeter observations sampling the Rayleigh-Jeans tail of the dust SED are crucial for constraining the dust absorption curve. During Early Science Cycle 2 (ES2) for the Large Millimeter Telescope (LMT), we carried out a 20 hour survey of the dust continuum at 1.1 mm on the central ≈ 200 pc of our Galaxy with the AzTEC bolometer array camera <cit.>. The AzTEC survey outperforms pre-existing FIR/submillimeter surveys (SPIRE/Herschel, Bolocam/CSO, HFI/Planck) 
  with regard to spatial resolution (HPBW=10.5”). Existing studies of the dust emission in the CMZ are mostly based on the Herschel Hi-GAL survey. Herschel/SPIRE (160-500 μ m) has comparatively low spatial resolution (HPBW_500μ m=36”) and insufficient spectral coverage of the Rayleigh-Jeans tail of the dust SED. Adding a high-resolution survey at 1.1 mm to the current data set significantly enhances our capability to uncover small scale structures and place tighter limits on beta <cit.>. 
  
  
  In this paper, we present the results of the AzTEC survey of the properties of dust in the CMZ, and study how parameters inferred from dust SEDs rely on different assumptions and priors. The structure of this paper is organized as follows. The observation strategy and data reduction for our AzTEC survey are briefly described in Section <ref>, also described in Section <ref> is how Herschel-SPIRE/Planck-HFI/CSO-Bolocam observations are processed and included into our study. A much more detailed discussion of both parts is presented in a separate paper on the data reduction of the AzTEC survey and the Bayesian analysis methodology of the multi-wavelength data (Tang et al.2020, hereafter Paper I). We extend the analysis to include the separation of the global background from the emission of the Galactic disk. In Section <ref> we describe our detailed SED analysis. The results are presented in Section <ref>. The scientific implications are discussed in Section <ref>. In Section <ref> we draw our conclusions.
  
  
  
  
  § OBSERVATION & DATA REDUCTION
  
  
  
  
   §.§ AzTEC 1.1 mm Survey of the CMZ
  
  
  The AzTEC survey was conducted during ES2 for the 32 meter LMT, from Apr 17 to June 18, 2014, with a total integration time of ≈ 20 hours. The survey covers the Galactic Center Region l=[-0.7, 0.9], b=[-0.6, 0.5], which roughly extends from Sgr B2 to Sgr C. The target field was mosaiced by square tiles, each observed with a raster-scan mode. The region immediately around Sgr B2 is excluded from this analysis due to Sgr B2's brightness, which exceeds 10-15 Jy and is subsequently not treated properly by our analysis due to the extreme dynamic range of flux and the possibility that the detector response is non-linear for the source. The optimal spatial resolution is 8.5” for the 32 meter LMT. However, since we adopted a high scanning speed of 200”/s, we eventually obtain a beam size of 10.5” as a result of under-sampling. The 1σ noise level is about 15mJy/beam after a 20 hours integration time.
  
  The raw data were reduced using the standard AzTEC analysis pipeline <cit.>. 
  We use iterative Principle Component Analysis (PCA) to remove correlated signals among bolometers, which are primarily contributed by the atmosphere, emissions from the telescope itself and non-Gaussian noises associated with the secondary mirror and back-end instruments.
  PCA Cleaning is performed iteratively, until a conversion is reached such that the rms in a final noise map is consistent with no significant astronomical signal.
  
  
  
  
  
  
  
  
  
   §.§ Processing of Herschel, Planck and CSO/Bolocam maps
  
  
  We create a combined 1.1 mm map from the AzTEC 1.1 mm map, the Planck/HFI 353 GHz map (Planck 2013 data release (PR1)) and the CSO/Bolocam 1.1 mm map <cit.>, to compensate for the large scale emission filtered out by the PCA cleaning in the AzTEC map. The Planck/HFI (353 GHz or 850 μ m) map is scaled to 1.1 mm to match the wavelength of the AzTEC & Bolocam maps before merging. We apply pixel-by-pixel colors corrections and scaling factors from 850 μ m to 1.1 mm to the Planck map inferred from best-fit 160 μ m-850 μ m SEDs. CO J=3-2 contamination has been removed from the Planck map, using estimates provided by the Planck Legacy Archive. The combined map is created using an approach described by <cit.>. This approach is mathematically equivalent to ”feathering”, an approach widely used for combining interferometer observations with single-dish observations, but is performed here in the “map domain” instead of the Fourier domain.
  
  The 1σ statistical noise is 15 mJy/beam in the AzTEC map, comparable to that in the Bolocam map. The noise in the Planck/HFI map is negligible. We further apply a 10% relative calibration uncertainty to the final 1.1 mm compound map, which accounts for beam variantions in the AzTEC maps. 
  
  To construct dust SEDs, we further take advantage of existing Herschel PACS/SPIRE 160 μ m, 250 μ m, 350 μ m and 500 μ m maps from the Hi-GAL survey <cit.>, which have been color-corrected using the Photometer Calibration Products from the ESA Herschel Science Archive. The errors in the Herschel maps are dominated by calibration uncertainties, which could be divided into relative calibration uncertainties and absolute calibration uncertainties. We adopt a relative calibration uncertainty of 2% for all SPIRE bands, and a relative uncertainty of 5% for the PACS 160 μ m band <cit.>. We notice that some authors adopted more conservative estimates of the relative uncertainties for extended sources, inferred from comparisons between observations taken by Herschel/PACS and those made by other facilities. (e.g. Spitzer/MIPS, AKARI, <cit.>). However, the fluctuations in the low surface brightness region of the PACS 160μ m maps indicate that the relative uncertainty should be <5%. Furthermore, it is dangerous to model absolute calibration offsets without an accurate knowledge of the dust absorption curve since these two ingredients are degenerate. Therefore, we ignore the absolute calibration uncertainties.
  
  
  
   §.§ Fore/Background Estimates
   
  
  
  
  
  
  
  In order to separate the column densities of the CMZ clouds, N_cmz, from the foreground and background column densities N_fb, we assume that along each line of sight, the dust emission from the fore/background is only a function of Galactic latitude and that, for each, the column density and temperature exponentially decrease away from the Galactic plane. 
  
  
      N_fb = N_fb,0× exp(- |b-b_0N|/σ_N )
  
  
  
      T_fb = T_fb,0× exp(- |b-b_0T|/σ_T )
  
  
  where N_fb,0 and T_fb,0 are peak column density and peak temperature, 
  b_0{N,T} and σ_{N,T} are offsets and
  scale heights, respectively. The spectral index β is fixed to 1.8 for fore/background dust emission, which is derived from SED-fitting to the Herschel 160-500 μ m and Planck 353 GHz maps degraded to the lowest resolution of the Planck map.
  
  The above model is constrained from three “pure” fore/background regions at high Galactic latitudes, which are shown in Figure <ref>. These low-flux regions are visually selected from the 500 μ m and 1.1 mm maps. N_fb,0, T_fb,0, b_0{N,T} and σ_{N,T} are derived by fitting 160-850 μ m dust SEDs from these three regions, assuming that there is no CMZ component, pixel by pixel, after degrading every map to the lowest resolution at 850 μ m. Figure <ref> shows a comparison between the fitted model and the observed flux densities in the three defined “pure” fore/background regions. The blue shaded area shows ± 1σ flucations of the observed intensities at different latitudes. The best-fit model has (N_fb,0[cm^-2]=10^22.08, b_0N=-0.05^∘, σ_N=0.66^∘) and (T_fb,0=24.8K, b_0T=-0.17^∘,σ_T=5.75^∘), notice that T_fb is almost constant across the region. σ_N=0.66^∘ corresponds to a scale height of 97 pc, consistent with previous findings <cit.>.
  
  This approach of the background subtraction is performed here in the {N,T} space rather than on each flux map, such as that used by <cit.>. In this way, we take advantage of the knowledge that flux densities in different bands are correlated to follow an approximated modified black-body SED. 
  
  
  
  § MODELING DUST PROPERTIES WITH BAYESIAN ANALYSIS
   
  
  In this section, we carry out a Bayesian analysis of the dust SEDs from the CMZ to explore the optical properties of dust grains. We adopt a forward modeling strategy to fit a dust model to multi-band maps, each diluted by a different instrumental PSF. The performance of this model-based deconvolution technique has been demonstrated in Paper I with a model of single temperature modified black-body (hereafter STMB). In this work, we extend and optimize this analysis to improve the estimation of physical parameters intrinsic to the CMZ.
  
  
  
   §.§ Single Temperature Dust Model
  
  Here we briefly recap the STMB model that is detailed in Paper I. This relies on three parameters: dust temperature T,  column density grid N_H_2 and dust spectral index β. The surface brightness F_i(ν_j) at pixel(i) and frequency ν_j is given by:
  
  
      F_i(ν_j) = [1-exp(-τ_i,ν_j)] B_ν_j(T_i) Ω_j
  
  
  where Ω_j is the beam area in the jth band. 
  B_ν_j(T_i) is the Planck function. τ_i,ν_j is the optical depth at frequency ν_j, which is given by:
  
  
      τ(i,ν_j) = κ_0 (ν_j/ν_0)^β_iμ m_H×N_H_2_i× 1%
  
  
  where κ_0 is the absorption cross section per unit mass at frequency ν_0. We adopt κ_0=1.37 cm^2/g and ν_0=c/1000 μ m from <cit.> for coagulated dust grains with thin ice mantles (their Table 1). We also adopt a mean molecular weight μ=2.8 per H_2 molecule from <cit.> and a dust-to-gas mass ratio of 1% to convert from N_H_2 to column dust mass density. This model is not restricted to an optically thin approximation (κ∝ν^β).
  
  The raw flux map 𝐅(ν_𝐣) calculated above is diluted to the instrumental resolution of each wavelength band to match the data:
  
      𝐌𝐨𝐝𝐞𝐥(ν_𝐣) = 𝐅(ν_𝐣)⊗𝐛𝐞𝐚𝐦_𝐣
  
  
  where ⊗ refers to convolution. All beams profiles are approximated as Gaussian profiles. The Full Width Half Maximum (FWHM) of the beams are 13.6” at 160 μ m, 23.4” at 250 μ m, 30.3” at 350 μ m, 42.5” at 500 μ m and 10.5” at 1.1 mm, respectively. The beam sizes of the PACS/SPIRE maps are larger than their nominal values <cit.>, which is due to the high scanning speed adopted by the Hi-Gal survey.
  
  
  
   §.§ STMB with Fore/Background Subtraction
  
  
  The total flux along any line of sight is the sum of the CMZ and its fore/background. With the self-absorption being accounted for, the total flux is:
  
  
      F_tot = F_bg× exp(-τ_cmz-τ_fg) + F_cmz× exp(-τ_fg) + F_fg
  
  
  where F_bg/F_cmz/F_fg are the intrinsic intensities of the background/CMZ/ foreground components, respectively. Because our concerned Galatic latitude range is small, we assume that the foreground and background intensities along each line of sight are identical: F_bg = F_fg = (1-exp(-1/2τ_fb)) B_ν(T_fb), τ_bg=τ_fg=1/2τ_fb.
  
  For MCMC sampling, it is more convenient to use the integrated column density, N_tot along each line of sight as a free parameter, instead of N_cmz. Then we have:
  
  
      0.5!F_tot
           = 
           F(N_tot, {T,β}_fb), if N_tot <= N_fb(a), 
      
      
           F({T,N,β}_bg) × exp(-τ_cmz-τ_fg) 
       
           +  F({T,N,β}_cmz) × exp(-τ_fg) 
      
           + F({T,N,β}_fg), if N_tot > N_fb(b),
  
  
  where F({T,N,β}) = [1-τ(N,β)] B_ν_i(T) while τ is the total internal opacity. N_cmz = N_tot - N_fb, τ_bg = τ_fg = 1/2τ_fb, N_bg = N_fg = 1/2N_fb, T_bg = T_fg = T_fb and β_bg = β_fg = 1.8. Note that fore/background fluxes are not completely fixed to the intensities calculated from Eq <ref> and Eq <ref>. In Eq <ref> (a), while T is always fixed to the values T_fb inferred from the best-fit Eq <ref> and β is to 1.8, N_tot could take values smaller than N_fb. In other words, low column density cells are not elevated to the values inferred from Eq <ref>.
  
  
  
   §.§ MCMC Analysis: Sampling Strategy
   
  
  In Paper I, we have demonstrated a Markov chain Monte Carlo (MCMC) approach for a model-based deconvolution of multi-band maps. A key ingredient of this procedure is a smoothness prior. In this section, we start by recapping the smoothness prior that is defined in Paper I, and then propose a more generalized form of smoothness prior that is applicable to the STMB model with the fore/background subtraction. 
  
  As is well known, deconvolution normally results in over-fitting to noise, which manifests as high-frequency fluctuations among neighboring cells in the best-fit maps. We employ regularized Bayesian inference <cit.> to mitigate this issue. Briefly, we adopt a simple form of smoothness prior based on the local gradients of the parameter x_ix,iy to be sampled:
  
  
      ln(P(x_ix,iy)) = ln(L(x_ix,iy)) + P_1
  
  
  in which:
  
  
      P_1 =  -1/2∑_j=-1, 1(x_ix+j,iy - x_ix,iy)^2 + ∑_k=-1, 1(x_ix,iy+k - x_ix,iy)^2/2σ_1^2
  
  
  Here, ln(P(x_ix,iy)) and ln(L(x_ix,iy)) are the logarithms of the full conditional posterior and likelihood for parameter x_ix,iy, where 𝐱 is {N_H_2,T,β} while P_1 is an a priori smoothness of 𝐱 and is estimated with gradients among all adjacent cells in a grid. σ_1 is user-defined and could be viewed as an a priori mean standard deviation of adjacent cell-cell differences. In Paper I, we explore the choices of σ_1 and show that, as demonstrated by spectral density analysis, for a broad range of σ_1, the performance of this regularized deconvolution approach is superior to that of the conventional approach that involves degrading every image to the lowest resolution. 
  
  For an STMB with the fore/background subtraction, however, the above prior has two issues. First, since there is no reason to assume any continuity between a CMZ component and a fore/background component, P_1 for T_cmz or β_cmz should be defined only among pairs of cells both having a CMZ component (i.e., N_tot>N_fb). And since P_1 has a negative contribution to the posterior, it yields a bias against having a CMZ component in every cell. Furthermore, in low density regions, where a cell having a CMZ component can have less than 4 neighbors which also have a CMZ component, Eq <ref> has less constraining power. To solve both issues, we propose a more generalized form of smoothness prior, calculated from average gradients upon a n × n block centered at each cell to be sampled:    
  
  
      P_n = -2 ∑_j=-n-1/2^n-1/2∑_k=-n-1/2^n-1/2 (x_ix+j, iy+k - x_ix,iy)^2 w_j,k/(2σ_n^2)/∑_j=-n-1/2^n-1/2∑_k=-n-1/2^n-1/2 f_cmz,j,k
  
  
  where n is an odd number, while w_j,k is defined as:
  
  
      w_j,k = 1/√(j^2+k^2), if N_tot,0,0>N_fb,0,0 and N_tot,j,k>N_fb,j,k, 
      
                0, if N_tot,0,0<=N_fb,0,0 or N_tot,j,k<=N_fb,j,k,
  
  
  and f_cmz is:
  
  
      f_cmz,j,k =  1, if N_tot,0,0>N_fb,0,0 and N_tot,j,k>N_fb,j,k, 
      
                0, if N_tot,0,0<=N_fb,0,0 or N_tot,j,k<=N_fb,j,k,
  
  
  Notice that from Eq <ref> to Eq <ref>, the scaling factor changes from 1/2 to 2. In Eq <ref>, the factor of 1/2 accounts for the fact that each pair is counted twice when P_1 is summed over all cells. In Eq <ref>, the factor of 2 ensures that σ_n is defined in a comparable fashion to σ_1, since gradients in Eq <ref> are summed over all adjacent pairs, which have a total number of ≈ 2N_cell.
  
  Throughout this study, we use both P_1 and P_n defined above to relieve over-fitting during forward modeling. P_n is adopted only for the STMB model with the fore/background subtraction. The smoothness priors for different models discussed in this paper are summarized in Table <ref>, where σ_1 corresponds to P_1 and σ_n (n>1) corresponds to P_n.
  
  
  
  
  
  
  
   §.§ Hierarchical Bayesian Model
   
  
  A common problem encountered in physical modeling is parameter degeneracy. The global distribution of the estimated parameters could be viewed as a convolution of their natural distribution with the probability distributions of their estimated values propagated from measurement uncertainties. Given that measurement uncertainty is always present, parameter degeneracy leads to correlated probability distributions, which dilute the apparent distribution of the best-fit parameters toward a false correlation. In SED analysis with a STMB model, T and β are known to have a high degeneracy, which manifests as a banana-shaped posterior distribution. This makes it difficult to recover the intrinsic T-β. <cit.> examined several existing techniques aiming to recover the intrinsic (T,β) relation and concluded that all techniques suffer from some bias. 
  
  Hierarchical Bayesian Analysis has been proposed to remedy the intrinsic correlation between T and β, by implementing the natural distribution of parameters as a prior to the model <cit.>.  Following <cit.>, we adopt a multivariate Student-t distribution as a prior for the T-β distribution. The posterior can be written as:
  
  
      0.45!P(μ,Σ | D) = ∏_i P(D | 𝐱_𝐢) P(𝐱_𝐢 | μ, Σ) × P(μ, Σ)
  
  
  
      0.5!P(𝐱_𝐢 | μ, Σ) ∝1/|Σ|^1/2× [1 + 1/d(𝐱_𝐢 - μ)^TΣ^-1 (𝐱_𝐢 - μ)]^-(d+2)/2
  
  
  
      𝐱_𝐢=(ln(T_i), β_i)
  
  
  where D is the data, and μ is the global mean of 𝐱_𝐢. Σ is the covariance matrix of 𝐱_𝐢. When d+1 samples are drawn from a normal distribution, Eq <ref> can be viewed as a distribution of the deviation of the sample mean from the true mean, divided by a normalized sample standard deviation. Here d is the degrees of freedom and is arbitrarily set to 8 following <cit.>. A Student-t distribution with smaller degrees of freedom has a larger portion of outliers relative to a normal distribution.
  
  The covariance matrix Σ can be decomposed as:
  
  
      Σ = SRS
  
  
  
  
  
  where S is the diagonal matrix of the standard deviations and R is the correlation matrix.
  
  With P(ln(T_i), β_i | μ, Σ) as an additional prior, we have 5 hyperparameters: 
  μ_ln(T), μ_β, σ_ln(T), σ_β and ρ_ln(T), β, which are sampled along with 
  (lg(N_i), ln(T_i), β_i) using a slice-within-Gibbs strategy. Since the covariance matrix Σ is a 2x2 matrix, 
  it is always positive-definite as long as -1< ρ_ln(T), β< 1. We can simply place
  a uniform prior on ρ_ln(T), β between -1 and 1. We further give uniform priors on the rest of parameters: 
  μ_ln(T[K])∈ (ln(5), ln(60)), μ_β∈ (0.5, 3.0), σ_ln(T[K])∈ (0.02, 0.4), and σ_β∈ (0.02, 0.4). 
  
  An illustration of this hierarchical Bayesian analysis is shown in Figure <ref>. Here, we simulate a sample of 22400 dust SEDs from a 3x3 multivariate normal distribution, with (μ_lg(N[cm^-2])=22.5, μ_ln(T[K])=ln(20), μ_β=2), (σ_lg(N[cm^-2])=0.2, σ_ln(T[K])=0.2, σ_β=0.1)[lg means base-10 log] 
  and (ρ_lg(N), ln(T)= -0.5, ρ_lg(N), β=0.5, ρ_ln(T), β=0.3. The signal-to-noise ratios are identical to our observed data. ln(T) and β are simulated to follow a positive correlation with ρ_ln(T), β=0.3. The best-fit distribution of ln(T) and β derived from a regular Bayesian analysis shows an apparent anti-correlation. For this simulation, we adopt a multivariate prior (2x2 multivariate normal distribution) that correctly characterizes the natural distribution (3x3 multivariate normal distribution) and is able to accurately recover the intrinsic correlation between T and β. We caution that the natural T-β distribution in molecular clouds might not follow a multivariate bell-shape distribution as we simulated. However, based on a magnetohydrodynamical (MHD) simulation of molecular clouds, <cit.> has shown that, multivariate prior could retain the information of the correlation coefficient under
  reasonable noise levels.
  
  
  
  § RESULTS: SINGLE TEMPERATURE MODIFIED BLACKBODY MODEL (STMB)
   
  
  
  
  
  
  The product maps of T, N_H_2 and β before and after fore/background subtraction are shown in Figure <ref>.
  Overall, the ranges of T and N_H_2
  are similar to those derived by <cit.> with DUSTEM. The temperatures of dense clouds are typically ≲ 20K, 
  and the peak column density N_H_2 is ≈ 10^23.5 cm^-2. Below N_H_2≈ 10^22 cm^-2, 
  the fluxes are dominated by fore/background emission. 
  
  As discussed in Paper I, the effective resolution in each best-fit map derived by our model-based deconvolution approach is between 10.5” and 42.5” and is different from parameter to parameter. While N_H_2 strongly depends on the combined 1.1 mm map, which has the highest resolution, T and β are more dependent on lower resolution maps. The final achieved resolution also relies on the choice of the smoothness prior. With a smaller cell size, the number of free parameters increases, but there is no extra information on sub-PSF scales. In other words, the effective degrees of freedom are reduced. Therefore, stronger smoothness priors are required to avoid overfitting. We have examined the best-fit maps after the fore/background subtraction with two different configurations, one configuration with a smaller cell/pixel size (7”) and strong priors (σ_1, lg(N)=0.05, σ_5,ln(T)=0.05, σ_5,β=0.05), a second one with a larger cell/pixel size (14”) and weak priors (σ_1, lg(N)=0.1, σ_5,β=0.1). A close-up comparison of the two results is shown in Figure <ref>. On scales larger than 42.5”, there is no apparent difference. Since this particular study focuses on the global distributions of dust properties in the CMZ, we determine to use a cell/pixel size of 14” for all models so as to avoid bias induced by the priors.
  
  On large scales, there are two pronounced correlations: a negative correlation between N_H_2 and T and a positive correlation between N_H_2 and β. 
  The marginalized distributions of the best-fit T, N_H_2 and β are plotted in Figure <ref>. In each panel we also plot three typical projected sampled posteriors at different locations in the parameter space. Cells in high galactic latitudes (b<-0.2^∘ or b>0.1^∘) are excluded. Measurement uncertainties are partially responsible for the apparent correlation between estimated T and β , which propagates into a banana-shaped posterior distribution. However, the sampled posterior distributions suggest that a genuine anti-correlation between T and β is present. The hierarchical Bayesian analysis also supports an intrinsic T-β anti-correlation. Figure <ref> shows the T-β distribution derived by modeling the natural T-β distribution as a multivariate Student-t prior distribution (Section <ref>). 
  The T-β distribution does not significantly change other than a reduction of high-temperature cells. 
  The estimated correlation coefficient ρ_ln(T), β=-0.74 indicates a strong anti-correlation. On the other hand, through the same analysis, we find that ρ_lg(N_H_2), β=0.89, indicating an even stronger correlation between N_H_2 and β.
  
  
  
  
  
  
  
  
  
  In Figure <ref> we plot the histograms of the ratios between the best-fit flux and the observed flux, F_fitted/F_observed. There is a systematical offset of ≈10% at 1.1 mm. This large systematic offset in contrast with Herschel bands is partially due to more substantial uncertainties at 1.1 mm. It is not clear whether this systematic offset is model-driven or due to calibration error.
  We do see a small systematic difference between high-density cells and low-density cells, which is however contrary to what is anticipated if the variation of β is related to the filtering effect in the Bolocam and the AzTEC maps at 1.1 mm. Indeed, dense clouds should be less affected by the filtering effect.
  As demonstrated in Figure <ref>, with the 1.1 mm map being removed, we perform the same analysis on the Herschel 160-500 μ m maps and find a distribution of β similar to but systematically smaller than that in Figure <ref>. 
  
  We also notice that some identified foreground objects in the CMZ show no sign of elevated β. <cit.> suggest that the dark cloud G359.94+0.17 is composed of two clouds in the foreground, 
  with V_lsr=0 km/s and 15 km/s. The comet-like feature near Sgr C complex (l=359.64, b=0.24)
  is associated with a foreground HII region RCW 137 <cit.> 1.8 kpc away.
  These two regions are marked in Figure <ref>. 
  
  
  
  The CMZ is moderately optically thick at 160μ m, partially due to high column densities in the CMZ and partially due to the steep slope of the dust absorption curve. The highest optical depth is τ_160≈ 1. τ_160 
  is irrelevant to our choice of the amplitude of κ_0, 
  since κ_0 and N_H_2 are completely degenerate. 
  
  
  
  § DISCUSSION
   
  
  
  
   §.§ Increased β in Dense Clumps
  
  
  
  
  We have identified a positive-correlation between N_H_2 and β. While this trend is qualitatively in agreement with existing observations <cit.>,  
  increased β up to 2.4 towards density peaks can not be easily
  explained by existing dust models. <cit.> reported such a steep absorption curve in the CMZ based on ISO observations. The origin of this trend deserves some discussion.
  
  We noticed that in a recent study of the dust in the CMZ, <cit.> reported a null detection of any correlation between T and β, based on an analysis of the 160-500μ m data from the Hi-Gal survey, same as what we use here. However, this null detection could be due to their coverage of a much larger region, extending from l=[-1.2, 1.2], b=[-0.5, 0.5], to both higher and lower latitudes and beyond SgrB2 in the positive Galactic longitude direction. Our study of correlations are confined to a smaller region of l=[-0.7, 0.6], b=[-0.2, 0.1]. In fact, their results show likely a noticeable positive correlation between N_H_2 and β (their Figure 5) in our region of interest.   
  
  Previous studies on the molecular cloud “Brick” <cit.> adopted a fixed β=1.2 for modeling dust emission, which is significantly lower than our results. 
  This low value of β was proposed based on a comparison between the Herschel 500 μ m map and the Atacama Large Millimeter/submillimeter Array (ALMA) 3 mm dust continuum by <cit.>, who find that by adopting β=1.2, the scaled Herschel 500 μ m map best recovers the missed large scale emission at 3 mm in the spatially filtered ALMA map. This comparison was not quantitatively detailed, and the uncertainty is not clear. Their use of lower β should be considered an assumption instead of a measurement. 
  Indeed, contrary to our results, a flattening of dust spectral index in the millimeter portion of the SED has been reported by  Herschel and Planck studies in some environments <cit.>. The origin of this flattening is not clearly understood; potential candidates are discussed in <cit.>, including 1) an extra cold dust component; 2) dust growth in very dense clouds; 3) magnetic dipole emission; and 4) low energy transitions in amorphous solids. 
  
  
  
  The total column densities we measured are similar to those derived from recent studies <cit.>. The major sources of bias/uncertainty in our analysis include the assumption of the metallicity, the single temperature approximation, and the variation of the dust absorption curve. The metallicity in the Galactic Center is probably twice higher than the solar metallicity that we have assumed here <cit.>. The column densities are likely underestimated by a factor of < 2 using a single temperature approximation (Tang et al.2020, in preparation). Finally, as we'll show later in Section 5.2, the variation of the dust absorption curve could potentially lead to an overestimate of the column densities by a factor of ≈ 2. 
  
  It is not a trivial task to recover the intrinsic T-β relation. <cit.> have discussed spurious correlation due to temperature mixing along the line of sight. However, this effect is more likely to suppress β with additional cold components on the Rayleigh-Jeans tail, which cannot explain the increase of β in the dense clouds. Our data sample the Rayleigh-Jeans tail down to 1.1 mm, where the spectral slope depends only weakly on the temperature. Therefore, we conclude that the observed anti-correlation is largely intrinsic, as confirmed by our hierarchical Bayesian analysis.
  
  
  
   §.§ Dust Model Predicting Higher β in Dense Regions
  
  
  The spectral index of dust absorption is expected to be environment dependent, e.g., due to  dust growth (via accretion & coagulation, <cit.>), or to dust destruction (e.g., shattering & sputtering, <cit.>). In dense molecular clouds, dust growth is usually expected due to high-frequency collision & sticking with low relative velocities. Classic models of dust growth suggest that this process leads to a lowering of β in submillimeter/millimeter wavelengths <cit.>. A recent model developed by <cit.> with updated optical properties of hydrogenated carbon grains could, however, reproduce the increase of β from FIR to submillimeter (≳ 500 μ m) by introducing a new population of small hydrogenated carbon grains. This results from a transition from aromatic-rich (i.e. hydrogen poor) hydrocarbons to aliphatic-rich (hydrogen-rich) ones in dense regions. For large grains, UV photo-process can at most aromatize down to a ≈ 20nm depth. Aliphatic-rich carbon grains have almost negligible emissivity in FIR-millimeter comparing to aromatic-rich ones and silicate grains. As a result, the spectral index of large grains is dominated by silicate features. Still, this model does not suggest β as high as ≳ 2. It is also questionable that dust growth could occur in the dense region in the CMZ, where the turbulent velocity dispersion is enhanced by a factor of a few <cit.>. Recently, <cit.> used DUSTEM to study the 3.6-70 μ m dust SEDs of the Arched Filaments in the CMZ, and suggest a depletion of large dust grains, which is in line with our finding that there is a millimeter deficit instead of an excess.
  
  
  
   §.§ Dust Model Predicting a β-T Anti-Correlation
  
  
  Laboratory experiments on “astrophysically relevant dust analogs” suggest complex relationships between the FIR-mm spectral index and the chemical composition or the physical structure (e.g., amorphous v.s. crystalline) of dust grains <cit.>. In these studies, an anti-correlation between T and β for amorphous dust is commonly reported. This correlation could also be reproduced by the TLS (two-level system) model proposed by <cit.>, who adopt a disordered charge distribution (DCD) on the nanometer scale and two-level systems on the atomic scale to describe the optical properties of dust. The absorption due to the DCD process is temperature independent and the combined absorption due to the TLS process, including resonant absorption, tunneling, and hopping, increases with temperature. This model is later applied by <cit.> to successfully reproduce the SEDs of ultracompact HII regions and cold clouds observed with Herschel/PACS & SPIRE and CSO/Bolocam. Both <cit.> and <cit.> reported an anti-correlation between T and β from large samples of cold clouds, which suggests that dust growth is at least not always a dominant factor in determining the spectral index. Our results confirm that this anti-correlation still exists in the more extreme CMZ environment.
  
  
  In the TLS model, the unnormalized absorption coefficient Q_abs can be divided into four components <cit.>:
  
      Q_abs = Q_DCD + A(Q_res+Q_phon+Q_hop)
   
  where A is a material-dependent parameter determining relative amplitudes of the temperature-independent DCD and the temperature-dependent TLS terms. The TLS terms, which become more important at long wavelengths, are further divided into three terms: resonant absorption Q_res, phonon-assisted tunneling relaxation Q_phon and hopping relaxation Q_hop. A simplified TLS model provided by <cit.> reduces the dust absorption curve to a function of only temperature T and wavelength λ: Q_abs=Q_abs(λ, T), with material-dependent parameters A, l_c, and c_Δ determined from FIR to millimeter SEDs, separately, for two samples: a sample representing the diffuse medium (FIRAS/WMAP) and a sample representing the cold dense environment (Archeops). The charge correlation length, l_c, controls asymptotic behaviors of Q_DCD, which approaches λ^-2 at short wavelengths and approaches λ^-4 at long wavelengths <cit.>. c_Δ is an additional parameter of tunneling states in Q_hop <cit.>. The values and uncertainties of A, l_c, and c_Δ in both environments are listed in Table 3 of  <cit.>. The corresponding absorption curves for the cold dense environment are shown in Figure <ref>. The power-indices β from 500-1100 μ m are in the range of 1.7-1.3 for the diffuse medium and 2.1-1.0 for the cold dense environment, decreasing with increasing temperature, from T=10-40 K.
  
  We can test to what extent the TLS model could be used to describe the dust emission in the CMZ. Here we focus on the region between -0.2^∘<b<0.1^∘, where the gas density is the highest, and adopt the TLS model for grain properties in the dense environments (solid lines in Figure <ref>). Figure <ref> shows the best-fit maps and illustrates the goodness of the fitting. The intensity deviation of the data from the best-fit model to the observations is most significant in the 1.1 mm band, with a factor close to 50% deficit. Apparently, the TLS model can not achieve higher β and lower temperature: β_500μ m-1.1mm is at most 2.1 at a temperature of 10K and 1.8 at 20K. At shorter wavelengths, between 100 - 300 μ m, the TLS model shows a very weak dependence on temperature. These deviations, which may partly due to the temperature
  mixing along the line of sight, might suggest a further difference between the grain properties in the CMZ and those in a typical Galactic dense environment. 
  
  
  
  
  The mass absorption coefficient of dust grains can be lowered at millimeter wavelengths if they have a crystalline structure <cit.>. In crystalline material, only a small number of phonons (lattice vibrations) can contribute to FIR absorption. The disorder of the atomic arrangement in amorphous materials leads to a breakdown of the selection rules for the frequency/wavenumber that govern the excitation of vibrational modes, which induces in the longest wavelength range a broad absorption band. This difference between crystalline and amorphous materials could be observed in a broad temperature range, between ≈ 10-300K <cit.>. However, the formation of crystalline dust grains usually requires condensation or annealing with T ≳ 1000 K, and crystalline silicate, with identifiable spectral features, and are thus expected to occur primarily in the circumstellar environment, occasionally in diffuse ISM where shocks are present <cit.>. Furthermore, it is expected that crystalline dust undergoes amorphization in the ISM environment <cit.>. In principle, strong shocks which prevail in the CMZ could potentially produce temperatures high enough for crystallization, but this scenario is yet to be explored.  
  We conclude that the observed T-β anti-correlation could be explained by four non-mutually-exclusive possibilities. 1) dust growth impediment and shatterings in the turbulent CGM environment; 2) an intrinsic T-β anti-correlation and 3) dust growth involved with hydrogenated carbons. A combination of 3) with either 1) or 2) remains the most plausible scenario.
  
  
  
  
  
  
  
  
   §.§ Impacts of β Variation
   
  
  Both the dust growth scenarios <cit.> and the TLS model <cit.> suggest a wavelength-dependent change of β. It is, therefore, worthwhile to explore the deviations from a single power-law absorption curve. Consider a smoothly broken power law for κ_λ:
  
  
      κ_λ = κ_λ_t (λ/λ_t)^-β_1{1/2[ 1 + (λ/λ_t)^1/δ] }^(β_1 -β_2)δ
  
  
  
  For simplicity, we fix δ to 0.1, this leads to a relatively sharp transition from β_1 to β_2 at wavelength λ_t. Limited by the five-bands of our SEDs, β_1 is almost completely degenerate with the dust temperature. We thus fix β_1 to a value of either 1.5 or 2.0 to investigate its impact on the measurement of the temperature and the column density. The term κ_0 (ν/ν_0)^β in Eq <ref> is replaced by Eq <ref>. The STMB model with this new β curve is fitted to the CMZ maps, again we apply smoothness priors to the parameter grids to avoid over-fitting. We set σ_lg(N_H_2)=0.1, σ_β_2=0.2 and σ_λ_t=20 μ m. The best-fit parameter distributions for the model with β_1=2.0 are shown in Figure <ref>. We find that β_2=2-3 and λ_t ≈ 500 μ m across the map, which is within the range of those suggested by recent experimental studies on astrophysically relevant dust analogs <cit.>, a summary of which is given in Table 1 of <cit.>. Figure <ref> provides a comparison between the temperature and column densities derived from the best-fit models with the single β and with a broken power-law β curve. This comparison illustrates that the estimation of the temperature is very sensitive to the assumption of κ_λ. The temperatures derived with β_1=1.5 are systematically higher by 20-50% compared to those derived with a β_1 =2.0, as a result of the reduced absorption coefficient at short wavelengths, which also leads to lower optical depths at 160 μ m. Still, this effect can not fully explain the discrepancy between the dust temperature and the gas temperature in the CMZ, which is a factor of 2-5 <cit.>.
  
  
  
  
  § SUMMARY
   
  
  To explore dust properties in the CMZ, we have combined the AzTEC 1.1 mm map with existing Herschel, Plank and Bolocam surveys from 160 μ m to 1.1 mm and carried out a joint SED analysis. We have developed an MCMC analysis tool which incorporates the knowledge of the PSFs to improve the spatial resolution, as well as the treatment of global background emission in different bands Equipped with this technique, we have explored the spatial variation of the column density, the dust temperature, and the dust spectral index β in the CMZ. Our main results and conclusions are the following:
  
  
  1) The spectral index β of the dust absorption curve increases from 1.8 to 2.4 from intermediate column densities (N_H_2≈ 22.5 cm^-2) to high densities (N_H2≈ 23.5 cm^-2). We confirm with a hierarchical Bayesian analysis that this correlation is not due to model degeneracy. We also derive a similar distribution of β by only using Herschel/Planck maps. Furthermore, we notice an absence of increased β toward foreground dense clouds in the same field. Therefore, the increase of β towards cold and dense clumps is induced by the CMZ environment.   
  
  
  2) The positive correlation between N_H_2 and β can be qualitatively, but not yet quantitatively, explained by contemporary dust models. This correlation could also be partially owing to a lack of dust growth, or even shattering due to the grain-grain collisions in a highly turbulent environment.
  In principle, the correlation could be caused by an intrinsic dependence of β on the temperature of dust. However, We find that the required dependence cannot be reproduced by either the dust growth model <cit.> or the TLS model <cit.>. 
  
  
  3) The inferred dust temperature is strongly dependent on the assumed dust absorption curve. We show that, different assumptions for β (λ < 500μ m) result in 0.1-0.2 dex difference in column and up to 50% difference
  in temperature. This model uncertainty is too small to be responsible for the decoupling between the gas temperature and the dust temperature observed in the CMZ (e.g., <cit.>). 
  
  
  
  § DATA AVAILABILITY
  
  The AzTEC data and products underlying this article are available at [<https://github.com/tangyping/products.git>]. The Herschel Hi-Gal products were provided by the Hi-Gal team by permission. The datasets from Planck telescope were derived from [<https://pla.esac.esa.int/>]. The datasets from CSO/Bolocam were derived from [<https://irsa.ipac.caltech.edu/data/BOLOCAM_GPS>].
  
  
  
  § ACKNOWLEDGEMENT
  
  
  We thank the referee for his constructive suggestions that improved this paper. The AzTEC instrument was built and operated through support from NSF grant 0504852 to the Five College Radio Astronomy Observatory. The authors gratefully acknowledge the many contributions of David Hughes in leading the LMT to its successful operational state. This work is partly supported by NASA via the grant NNX17AL67G. This work is also based on observations made with Herschel, Planck and the Caltech Submillimeter Observatory (CSO) telescope. Herschel is a European Space Agency cornerstone mission with science instruments provided by European-led Principal Investigator consortia and with significant participation by NASA. Planck is a project of the European Space Agency with instruments funded by ESA member states, and with special contributions from Denmark and NASA. CSO was operated by the California Institute of Technology under cooperative agreement with the National Science Foundatio.
  
  
  
  [Agladze et al.(1996)]agladze96 Agladze, N. I., Sievers, A. J., Jones, S. A., Burlitch, J. M., Beckwith, S. V. W., et al., 1996, ApJ, 462, 1026
  
  [Aguirre et al.(2011)]aguirre11 Aguirre, J. E., Ginsburg, A. G., Dunham, M. K., et al. 2011, ApJS, 192, 
  4 
  
  [Aharonian et al.(2006)]aharonian06 Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, , 439, 695
  
  [Arendt et al.(2019)]arendt19 Arendt, R. G., Staguhn, J., Dwek, E., Morris, M. R., Yusef-Zadeh, F., Benford, D. J., Kovács, A., Gonzalez-Quiles, J., arXiv:1909.08025
  
  [Bally et al.(1987)]bally87 Bally, J., Stark, A. A., Wilson, R. W., & Henkel, C. 1987, , 65, 13 
  
  [Balog et al.(2014)]balog14 Balog, Z., Müller, T., Nielbock, M., et al. 2014, Experimental Astronomy, 37, 129 
  
  [Battersby et al.(2011)]battersby11 Battersby, C., Bally, J., Ginsburg, A., et al. 2011, , 535, A128 
  
  [Bendo et al.(2013)]bendo13 Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, , 433, 3062 
  
  [Blain et al.(2002)]blain02 Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, , 369, 111 
  
  
  [Boudet et al.(2005)]boudet05 Boudet, N., Mutschke, H., Nayral, C., et al. 2005, , 633, 272 
  
  [Casey(2012)]casey12 Casey, C. M. 2012, , 425, 3094 
  
  [Casey et al.(2014)]casey14 Casey, C. M., Narayanan, D., & Cooray, A. 2014, , 541, 45 
  
  [Chapin et al.(2011)]chapin11 Chapin, E. L., Chapman, S. C., Coppin, K. E., et al. 2011, , 411, 505 
  
  [Chen et al.(2016)]chen16 Chen, M. C.-Y., Di Francesco, J., Johnstone, D., et al. 2016, , 826, 95
  
  [Clark et al.(2013)]clark13 Clark, P. C., Glover, S. C. O., Ragan, S. E., Shetty, R., & Klessen, R. S. 2013, , 768, L34 
  
  [Coupeaud et al.(2011)]coupeaud11 Coupeaud, A., Demyk, K., Meny, C., et al. 2011, , 535, A124 
  
  [Crocker et al.(2010)]crocker10 Crocker, R. M., Jones, D. I., Melia, F., Ott, J., & Protheroe, R. J. 2010, , 463, 65 
  
  [Crocker & Aharonian(2011)]crocker11 Crocker, R. M., & Aharonian, F. 2011, Physical Review Letters, 106, 101102 
  
  [Deguchi et al.(2012)]deguchi12 Deguchi, S., Tafoya, D., & Nagisa, S. 2012, , 64, 28 
  
  [Demyk et al.(2017A)]demyk17A Demyk, K., Meny, C., Lu, X.-H., et al. 2017, , 600, A123 
  
  [Demyk et al.(2017B)]demyk17B Demyk, K., Meny, C., Leroux, H., et al. 2017, , 606, A50 
  
  [Demyk et al.(2013)]demyk13 Demyk, K., Meny, C., Leroux, H., et al. 2013, Proceedings of The Life Cycle of Dust in the Universe: Observations, Theory, and Laboratory Experiments (LCDU2013). 18-22 November, 2013. Taipei, Taiwan.
  
  [Draine & Salpeter(1979)]draine79 Draine, B. T., & Salpeter, E. E. 1979, , 231, 438 
  
  [Draine(2006)]draine06 Draine, B. T. 2006, , 636, 1114 
  
  [Dupac et al.(2003)]dupac03 Dupac, X., Bernard, J.-P., Boudet, N., et al. 2003, , 404, L11 
  
  [Faridani et al.(2018)]faridani18Faridani, S., Bigiel, F., Flöer, L., Kerp, J., Stanimirović, S., 2018, AN, 339, 87
  
  [Federrath et al.(2016)]federrath16 Federrath, C., Rathborne, J. M., Longmore, S. N., et al. 2016, , 832, 143 
  
  [Galliano(2018)]galliano18 Galliano, F. 2018, , 476, 1445 
  
  [Ginsburg et al.(2013)]ginsburg13 Ginsburg, A., Glenn, J., Rosolowsky, E., et al. 2013, , 208, 14 
  
  [Ginsburg et al.(2016)]ginsburg16 Ginsburg, A., Henkel, C., Ao, Y., et al. 2016, , 586, A50 
  
  [Goldsmith et al.(1997)]goldsmith97 Goldsmith, P. F., Bergin, E. A., & Lis, D. C. 1997, , 491, 615 
  
  [Gordon et al.(2014)]gordon14 Gordon, K. D., Roman-Duval, J., Bot, C., et al. 2014, , 797, 85 
  
  [Hankins et al.(2017)]hankins17 Hankins, M. J., Lau, R. M., Morris, M. R., & Herter, T. L. 2017, , 837, 79 
  
  [Henning & Mutschke (1997)]henning97 Henning T., & Mutschke H., 1997, A&A, 327, 743
  
  [Heyer et al.(2018)]heyer18 Heyer, M., Wilson, G. W., Gutermuth, R., et al. 2018, MNRAS, 473, 2222 
  
  [Hirashita & Yan(2009)]hirashita09 Hirashita, H., & Yan, H. 2009, , 394, 1061 
  
  [Indriolo et al.(2014)]indriolo14 Indriolo, Nick, Neufeld, D. A., Gerin, M., PRISMAS, WISH , 2015, , 800, 40 
  
  [Jones et al.(2011)]jones11 Jones, D. O., West, A. A., Foster, J. B., et al., 2011, AJ, 142, 44
  
  [Jones et al.(2013)]jones13 Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, , 558, A62 
  
  [Juvela et al.(2013)]juvela13 Juvela, M., Montillaud, J., Ysard, N., & Lunttila, T. 2013, , 556, A63 
  
  [Juvela et al.(2015)]juvela15 Juvela, M., Demyk, K., Doi, Y., et al. 2015, , 584, A94 
  
  [Kauffmann et al.(2008)]kauffmann08 Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., II, Lee, C. W., 2008, A&A, 487, 993
  
  [Kauffmann et al.(2017A)]kauffmann17 Kauffmann, J., Pillai, T., Zhang, Q., et al. 2017, , 603, A89 
  
  [Kelly et al.(2012)]kelly12 Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, , 752, 55 
  
  [Kemper et al.(2004)]kemper04 Kemper, F., Vriend, W. J., & Tielens, A. G. G. M. 2004, , 609, 826 
  
  [Köhler et al.(2012)]kohler12 Köhler, M., Stepnik, B., Jones, A. P., et al. 2012, , 548, A61 
  
  [Köhler et al.(2015)]kohler15 Köhler, M., Ysard, N., & Jones, A. P. 2015, , 579, A15
  
  [Krieger et al.(2017)]krieger17 Krieger, N., Ott, J., Beuther, H., et al. 2017, , 850, 77 
  
  [Kruegel & Siebenmorgen(1994)]kruegel94 Kruegel, E., & Siebenmorgen, R. 1994, , 288, 929 
  
  [Kruijssen & Longmore (2013)]kruijssen13 Kruijssen, J. M. D., Longmore, S. N., 2013, , 435, 2598
  
  [Kwon et al.(2009)]kwon09 Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, , 696, 841 
  
  [Li et al.(2018)]li18 Li, L., Shen, S., Hou, J., Yuan, H., Xiang, M., et al, 2018, ApJ, 858, 75
  
  [Lis & Menten(1998)]lis98 Lis, D. C., & Menten, K. M. 1998, , 507, 794 
  
  [Longmore et al.(2012)]longmore12 Longmore, S. N., Rathborne, J., Bastian, N., et al. 2012, , 746, 117 
  
  [Magnelli et al.(2012)]magnelli12 Magnelli, B., Lutz, D., Santini, P., et al. 2012, , 539, A155 
  
  [Marsh et al.(2016)]marsh16 Marsh, K. A., Ragan, S. E., Whitworth, A. P., & Clark, P. C. 2016, , 461, L16 
  
  [Mennella et al.(1998)]mennella98Mennella, V., Brucato, J. R., Colangeli, L., Palumbo, P., Rotundi, A., Bussoletti, E., et al., 1998, ApJ, 496, 1058
  
  [Meny et al.(2007)]meny07 Meny, C., Gromov, V., Boudet, N., et al. 2007, , 468, 171 
  
  [Mills(2017)]mills17 Mills, E. A. C. 2017, arXiv:1705.05332 
  
  [Morris & Serabyn(1996)]morris96 Morris, M., & Serabyn, E. 1996, ARA&A, 34, 645 
  
  [Morris(2015)]morris15 Morris, M. R. 2015, Lessons from the Local Group: A Conference in honor of David Block and Bruce Elmegreen, 391 
  
  [Molinari et al.(2010)]molinari10 Molinari, S., Swinyard, B., Bally, J., et al. 2010, , 122, 314 
  
  [Molinari et al.(2011)]molinari11 Molinari, S., Bally, J., Noriega-Crespo, A., et al. 2011, , 735, L33 
  
  [Najarro et al.(2009)]najarro09 Najarro, F., Figer, D. F., Hillier, D. J., Geballe, T. R., & Kudritzki, R. P. 2009, , 691, 1816 
  
  [Neal R. M. (2003)]neal03 Neal, R. M. 2003, The Annals of Statistics, 31, 705
  
  [Odegard et al.(2016)]odegard16 Odegard, N., Kogut, A., Chuss, D. T., et al. 2016, , 828, 16
  
  [Oka et al.(2007)]oka07 Oka, T., Nagai, M., Kamegai, K., Tanaka, K., & Kuboi, N. 2007, , 59, 15 
  
  [Oka et al.(2019)]oka19 Oka, T., Geballe, T. R., Goto, M., Usuda, T., Benjamin, McCall, J., Indriolo, N., , 883, 54
  
  [Ossenkopf & Henning(1994)]ossenkopf94 Ossenkopf, V., & Henning, T. 1994, , 291, 943 
  
  [Paradis et al.(2011)]paradis11 Paradis, D., Bernard, J.-P., Mény, C., & Gromov, V. 2011, , 534, A118 
  
  [Paradis et al.(2014)]paradis14 Paradis, D., Mény, C., Noriega-Crespo, A., et al. 2014, , 572, A37 
  
  [Paglione et al.(1998)]paglione98 Paglione, T. A. D., Jackson, J. M., Bolatto, A. D., & Heyer, M. H. 1998, , 493, 680 
  
  [Pillai et al.(2015)]pillai15 Pillai, T., Kauffmann, J., Tan, J. C., Goldsmith, P. F., Carey, S. J., Menten, K. M., 2015, , 799, 74
  
  [Planck Collaboration et al.(2011)]planck11 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, , 536, A17 
  
  [Popping et al.(2017)]popping17 Popping, G., Somerville, R. S., & Galametz, M. 2017, , 471, 3152 
  
  [Rathborne et al.(2014)]rathborne14 Rathborne, J. M., Longmore, S. N., Jackson, J. M., et al. 2014, , 795, L25 
  
  [Rathborne et al.(2015)]rathborne15 Rathborne, J. M., Longmore, S. N., Jackson, J. M., et al. 2015, , 802, 125 
  
  [Rowlands et al.(2014)]rowlands14 Rowlands, K., Gomez, H. L., Dunne, L., et al. 2014, , 441, 1040 
  
  [Russeil et al.(2003)]russeil03 Russeil, D., Juvela, M., Lehtinen, K., Mattila, K., & Paatero, P. 2003, , 409, 135 
  
  [Sadavoy et al.(2016)]sadavoy16 Sadavoy, S. I., Stutz, A. M., Schnee, S., et al. 2016, , 588, A30 
  
  [Santini et al.(2010)]santini2010 Santini, P., Maiolino, R., Magnelli, B., et al. 2010, , 518, L154 
  
  [Schnee et al.(2014)]schnee14 Schnee, S., Mason, B., Di Francesco, J., et al. 2014, MNRAS, 444, 2303 
  
  [Scott et al.(2008)]scott08 Scott, K. S., Austermann, J. E., Perera, T. A., et al. 2008, , 385, 2225 
  
  [Shetty et al.(2009)]shetty09 Shetty, R., Kauffmann, J., Schnee, S., Goodman, A. A., & Ercolano, B. 2009, , 696, 2234 
  
  [Shetty et al.(2012)]shetty12 Shetty, R., Beaumont, C. N., Burton, M. G., Kelly, B. C., & Klessen, R. S. 2012, , 425, 720 
  
  [Shields & Ferland(1994)]shields94 Shields, J. C., & Ferland, G. J. 1994, , 430, 236 
  
  [Tanaka et al.(2014)]tanaka14 Tanaka, K., Oka, T., Matsumura, S., Nagai, M., & Kamegai, K. 2014, , 783, 62 
  
  [Traficante et al.(2011)]traficante11 Traficante, A., Calzoletti, L., Veneziani, M., Ali, B., de Gasperis, G., et al., 2011, MNRAS, 416, 2932
  
  [Warren & Dye(2003)]warren03 Warren, S. J., & Dye, S. 2003, ApJ, 590, 673 
  
  [Wilson et al.(2008)]wilson08 Wilson, G. W., Austermann, J. E., Perera, T. A., et al. 2008, , 386, 807 
  
  [Wright et al.(2016)]wright16 Wright, C. M., Do Duy, T., & Lawson, W. 2016, , 457, 1593 
  
  [Ysard et al.(2012)]ysard12 Ysard, N., Juvela, M., Demyk, K., et al. 2012, , 542, A21 
  
  [Ysard et al.(2013)]ysard13 Ysard, N., Abergel, A., Ristorcelli, I., et al. 2013, , 559, A133 
  
  
  
  
  