# The evolution of temperature and density structures of OB cluster-forming molecular clumps

Y. Lin, F. Wyrowski, H. Liu, A. Izquierdo, T. Csengeri, S. Leurini, K. Menten

► To cite this version:

Y. Lin, F. Wyrowski, H. Liu, A. Izquierdo, T. Csengeri, et al.. The evolution of temperature and density structures of OB cluster-forming molecular clumps. Astronomy and Astrophysics - A&A, 2022, 658, pp.A128. 10.1051/0004-6361/202142023. hal-03799850

HAL Id: hal-03799850

https://hal.science/hal-03799850

Submitted on 6 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The evolution of temperature and density structures of OB cluster-forming molecular clumps

Y. Lin<fnref n="1" />,*, F. Wyrowski<fnref n="1" />, H. B. Liu<fnref n="2" />, A. F. Izquierdo<fnref n="3" />,<fnref n="4" />, T. Csengeri<fnref n="5" />, S. Leurini<fnref n="6" />, and K. M. Menten<fnref n="1" />

Received 15 August 2021 / Accepted 2 December 2021

## ABSTRACT

**Context.** OB star clusters originate from parsec-scale massive molecular clumps, while individual stars may form in $\le 0.1$ pc scale dense cores. The thermal properties of the clump gas are key factors governing the fragmentation process, and are closely affected by gas dynamics and feedback of forming stars.

**Aims.** We aim to understand the evolution of temperature and density structures on the intermediate-scale ($\le 0.1–1$ pc) extended gas of massive clumps. This gas mass reservoir is critical for the formation of OB clusters, due to their extended inflow activities and intense thermal feedback during and after formation.

**Methods.** We performed ~0.1 pc resolution observations of multiple molecular line tracers (e.g., CH$_3$CCH, H$_2$CS, CH$_3$CN, CH$_3$OH) that cover a wide range of excitation conditions, toward a sample of eight massive clumps. The sample covers different stages of evolution, and includes infrared-weak clumps and sources that are already hosting an HII region, spanning a wide luminosity-to-mass ratio ($L/M$) range from ~1 to ~100 ($L_\odot/M_\odot$). Based on various radiative transfer models, we constrain the gas temperature and density structures and establish an evolutionary picture, aided by a spatially dependent virial analysis and abundance ratios of multiple species.

**Results.** We determine temperature profiles varying in the range 30–200 K over a continuous scale, from the center of the clumps out to 0.3–0.4 pc radii. The clumps' radial gas density profiles, described by radial power laws with slopes between −0.6 and ~−1.5, are steeper for more evolved sources, as suggested by results based on dust continuum, representing the bulk of the gas (~$10^4$ cm$^{-3}$), and on CH$_3$OH lines probing the dense gas ($\gtrsim 10^6–10^8$ cm$^{-3}$) regime. The density contrast between the dense gas and the bulk gas increases with evolution, and may be indicative of spatially and temporally varying star formation efficiencies. The radial profiles of the virial parameter show a global variation toward a sub-virial state as the clump evolves. The linewidths probed by multiple tracers decline with increasing radius around the central core region and increase in the outer envelope, with a slope shallower than the case of the supersonic turbulence ($\sigma_v \propto r^{0.5}$) and the subsonic Kolmogorov scaling ($\sigma_v \propto r^{0.33}$). In the context of evolutionary indicators for massive clumps, we also find that the abundance ratios of [CCH]/[CH$_3$OH] and [CH$_3$CN]/[CH$_3$OH] show correlations with clump $L/M$.

**Key words.** ISM: clouds – dust, extinction – ISM: structure – ISM: abundances – ISM: molecules – ISM: kinematics and dynamics

## 1. Introduction

Massive star-forming clumps are progenitors of OB clusters (Williams et al. 2000; Motte et al. 2018). They typically have masses of $\gtrsim 1000 M_\odot$ over a spatial scale of ~1 pc. Fragmentation and accretion processes of OB star clusters are strongly influenced by the evolution of the kinematics, density, and temperature structure of parsec-scale gas clumps (Girichidis et al. 2011; Lee & Hennebelle 2019; Padoan et al. 2020), and vice versa (Krumholz et al. 2012; Offner et al. 2009; Hennebelle et al. 2020). In particular, the stellar initial mass function (IMF) appears to be universal; it varies weakly from one environment to another in the Milky Way, indicating that the formation of

the most massive stars is deterministic, favoring particular physical environments instead of randomly occurring in molecular clouds (Kroupa et al. 2013). This, together with dominant feedback caused by massive stars, may determine the evolutionary track of massive clumps. Accordingly, observational evidence can be collected by sampling a wide range of clumps at different evolutionary stages.

The process of gas collapse resulting in protostars has been studied for decades; among the first are the works by Larson (1969), Penston (1969), and Shu (1977). These are commonly referred to as “outside-in” and “inside-out” models, indicating the succession of the spherical collapse of isothermal clouds, which describe the gas flows immediately prior to and after development of a protostar singularity, respectively. The density profiles of Larson (1969) and Penston (1969) exhibit a $r^{-2}$ relation, while the density profile of the inner free-falling and

<footnote n="1">¹ Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany
e-mail: ylin@mpifr-bonn.mpg.de</footnote>

<footnote n="2">² Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, ROC</footnote>

<footnote n="3">³ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany</footnote>

<footnote n="4">⁴ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="5">⁵ OASU/LAB-UMR5804, CNRS, Université Bordeaux, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="6">⁶ INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy</footnote>

<footnote>* Member of the International Max-Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.</footnote>
outer static envelopes of the Shu (1977) model follow $r^{-1.5}$ and $r^{-2}$, separately. On the other hand, when turbulent pressure is taken into account to explain the observed linewidth-size scaling relation, the logatropic (nonisothermal) gas follows a flatter profile proportional to $r^{-1}$ (e.g., McLaughlin & Pudritz 1997) in the outer region. Recently, the process of spherically symmetric cloud collapse has been revisited extensively: Coughlin (2017) presents solutions for arbitrary initial density profiles, extending to non-self-similar regimes, and Murray & Chang (2015) and Murray et al. (2017) consider turbulent energy generated from gravitational collapse and a highly dynamic system (compared to hydrostatic equilibrium assumed by Shu 1977). Furthermore, due to significant heating sources and high opacities associated with massive star-forming clouds, the assumption of isothermality might break down and a polytropic equation of state (EOS) needs to be introduced which quantifies the balance of gas cooling and heating and can incorporate turbulent pressure (Curry & McKee 2000). The polytropic index $\gamma$ (with $T \propto \rho^{\gamma-1}$) has been shown to have a decisive effect on the dynamical evolution of molecular clouds (Passot & Vázquez-Semadeni 1998; Spaans & Silk 2000) and eventually on the IMF (e.g., Klessen et al. 2007; Jappsen et al. 2005; Lee & Hennebelle 2018; Hennebelle et al. 2020). However, a simple analytical form of gas EOS may not fully capture the gas thermal properties when stellar radiation feedback sets in. The radiative feedback has been suggested to have an impact in determining the characteristics of the stellar spectrum (e.g. Bate 2009; Urban et al. 2010; Guszejnov et al. 2016, 2017). Such impact can be only moderate in setting the peak of the IMF in typical galactic environment (Hennebelle et al. 2020), but it may be critical for boosting the formation of the more massive stars (Krumholz et al. 2007; Hennebelle et al. 2020). Given these theoretical developments, the time is right to measure with observations the detailed gas temperature and density distribution inside massive clumps.

Understanding how the mass of massive clumps is distributed over different density regimes is fundamental to understanding the evolution of the star formation rate (SFR) and star formation efficiency (SFE) on larger physical scales (e.g., Lee et al. 2015; Parmentier 2019). On cloud scales ($\gtrsim$10 pc), the gas column density distribution follows a log-normal probability function (N-PDF) in a turbulent medium, while it develops a power-law tail as significant gravitational collapse commences in high-density regimes (e.g., Klessen 2000; Kritsuk et al. 2011). The relevant scales are readily resolved in nearby star-forming clouds and OB cluster-forming regions (e.g., Kainulainen et al. 2009; Lin et al. 2016, 2017). The power-law shape is thought to originate from power-law density profiles (Federrath & Klessen 2013, Myers 2015). Hence, measurements of clump density profiles can provide insights into how the dense gas of molecular clouds leads to the power-law excess of N-PDFs.

Most previous works on the density structure of massive clumps are based on single-dish observations of both continuum and molecular lines. Works discussing samples of sources include (but are not limited to) van der Tak et al. (2000), Mueller et al. (2002), Beuther et al. (2002), Hatchell & van der Tak 2003, Rolffs et al. (2011), Williams et al. (2005), and Palau et al. (2014). With the advent of wide-band receivers, especially those equipping interferometers, spatially resolved multi-line observations have become efficient (e.g., Beuther et al. 2007, Li et al. 2019, Gieser et al. 2021), and are indispensable to measure the broad density and temperature ranges associated with massive star formation.

We conducted a pilot survey of eight massive clumps with the Submillimeter Array (SMA) and the APEX telescope. For

the target clump selection we followed the criterion elaborated in Sect. 2.1. The main molecular lines of interest are listed in Table 1, which include multiple efficient thermometers and densitometers for massive clumps, as suggested by single-dish observations toward a statistically large sample (Giannetti et al. 2017; Leurini et al. 2004, 2007). We use various modeling methods to quantify the clump density and temperature structure using these lines. Throughout the paper, we follow the existing nomenclature in the literature (e.g., Williams et al. 2000; Zhang et al. 2009; Liu et al. 2012; Motte et al. 2018). In this way, massive molecular clumps refer to structures with sizes of ~0.5–1 pc, massive molecular cores refer to the <0.1 pc size structures embedded within a clump, and condensations refer to the distinct molecular substructures within a core. In Fig. 1, we provide a schematic picture of different scales of a molecular cloud. The physical characteristics across the scales, as elaborated above, are indicated for individual structures. We are interested in understanding the physical structure of massive clumps, which have a vast range of gas densities and are the building blocks of the star-forming clouds; in particular, they compose the high-density tail of the cloud N-PDF.

The paper is organized as follows. In Sect. 2, we describe the observations and data reduction. In Sect. 3, we describe the radiative transfer modeling procedure and elaborate on the analysis of both line and continuum observations to derive the temperature, density, linewidth, and virial parameter profiles, and the abundances of multiple species. Section 3.1 gives a general outline of the radiative modeling methods and procedures we adopted. Section 3.3 provides an overview of the distribution of the molecular lines used as thermometers and densitometers in this paper. Sections 3.4 and 3.5, in addition to Appendices E and F, focus on the radiative transfer modeling procedures and results of continuum and molecular lines. In Sect. 3.2, the properties of the sources extracted from SMA 1.2 mm continuum are presented. An analysis of some complementary lines is presented in Appendix D. In Sect. 4, we discuss the outcome of the modeling results, with a comparison between target sources and the physical implications. Finally in Sect. 5, we conclude on our findings.

## 2. Observations and data reduction

### 2.1. Source selection

The target sources are selected from the APEX Telescope Large Survey of the Galaxy (ATLASGAL; Schuller et al. 2009), and are located at a distance of 4–6 kpc (Urquhart et al. 2018). They cover different evolutionary stages, suggested by different $L/M$ (Fig. 2) and different signposts of star formation activity (see further details in Appendix A). For comparison, in Fig. 2 we present the distribution of luminosity and mass for ATLASGAL sources that (1) are located at a distance within 4–8 kpc with a radius of <2 pc (Urquhart et al. 2018), and (2) have masses and peak fluxes higher than 300 $M_\odot$ and peak flux $\sim$2 Jy beam$^{-1}$, respectively. The background contours illustrate the distribution of ATLASGAL sources in a distance range of 4–8 kpc, with masses over 300 $M_\odot$ and peak flux $\gtrsim$2 Jy beam$^{-1}$ with a radius of less than 2 pc.

In Fig. 2, we also include several evolutionary tracks. The empirical evolutionary tracks for massive clumps with envelope masses of 80, 140, 350, 700, and 2000 $M_\odot$ which will form a single protostar with varying accretion rates based on a turbulent core model (McKee & Tan 2003) are shown (gray lines), as derived in Molinari et al. (2008). Evolutionary tracks of massive
<table><caption>**Table 1.** Molecular lines of interest covered by the SMA observations.</caption><thead><tr><th>Transitions</th><th>Rest frequency<br>(GHz)</th><th><i>E</i><sub>up</sub><br>(K)</th><th>Transitions</th><th>Rest frequency<br>(GHz)</th><th><i>E</i><sub>up</sub><br>(K)</th></tr></thead><tbody><tr><td colspan="6" style="text-align:center; font-weight:bold;">(RxA lower sideband)</td></tr><tr><td>C<sup>34</sup>S 5–4</td><td>241.016</td><td>27.8</td><td>CH<sub>3</sub>OH <i>v</i><sub>t</sub> = 1 5<sub>-1,5</sub>-4<sub>-1,4</sub><i>E</i></td><td>241.203</td><td>326.2</td></tr><tr><td>CS 5–4</td><td>244.935</td><td>35.3</td><td>CH<sub>3</sub>OH <i>v</i><sub>t</sub> = 1 5<sub>2,4</sub>-4<sub>2,3</sub><i>A</i></td><td>241.192</td><td>333.4</td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">(RxA upper sideband)</td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>0</sub>–11<sub>0</sub></td><td>205.081</td><td>64.0</td><td>CH<sub>3</sub>OH <i>v</i><sub>t</sub> = 1 5<sub>0,5</sub>-4<sub>0,4</sub><i>E</i></td><td>241.206</td><td>333.5</td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>1</sub>–11<sub>1</sub></td><td>205.076</td><td>71.1</td><td>CH<sub>3</sub>OH <i>v</i><sub>t</sub> = 1 5<sub>3,3</sub>-4<sub>3,2</sub><i>E</i></td><td>241.180</td><td>357.4</td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>2</sub>–11<sub>2</sub></td><td>205.065</td><td>92.5</td><td>CH<sub>3</sub>OH <i>v</i><sub>t</sub> = 1 5<sub>-3,2</sub>-4<sub>-3,1</sub><i>E</i></td><td>241.180</td><td>357.4</td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>3</sub>–11<sub>3</sub></td><td>205.045</td><td>128.2</td><td>CH<sub>3</sub>CN 13<sub>0</sub>–12<sub>0</sub></td><td>239.138</td><td>80.3</td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>4</sub>–11<sub>4</sub></td><td>205.018</td><td>178.2</td><td>CH<sub>3</sub>CN 13<sub>1</sub>–12<sub>1</sub></td><td>239.133</td><td>87.5</td></tr><tr><td>H<sub>2</sub>CS 6<sub>0,6</sub>–5<sub>0,5</sub></td><td>205.987</td><td>34.6</td><td>CH<sub>3</sub>CN 13<sub>2</sub>–12<sub>2</sub></td><td>239.119</td><td>108.9</td></tr><tr><td>H<sub>2</sub>CS 6<sub>2,5</sub>–5<sub>2,4</sub></td><td>206.054</td><td>87.3</td><td>CH<sub>3</sub>CN 13<sub>3</sub>–12<sub>3</sub></td><td>239.096</td><td>144.6</td></tr><tr><td>H<sub>2</sub>CS 6<sub>2,4</sub>–5<sub>2,3</sub></td><td>206.159</td><td>87.3</td><td>CH<sub>3</sub>CN 13<sub>4</sub>–12<sub>4</sub></td><td>239.064</td><td>194.6</td></tr><tr><td>H<sub>2</sub>CS 6<sub>3,4</sub>–5<sub>3,3</sub></td><td>206.053</td><td>153.1</td><td>CH<sub>3</sub>CN 13<sub>5</sub>–12<sub>5</sub></td><td>239.023</td><td>258.9</td></tr><tr><td>H<sub>2</sub>CS 6<sub>3,3</sub>–5<sub>3,2</sub></td><td>206.053</td><td>153.1</td><td>CH<sub>3</sub>CCH 14<sub>0</sub>–13<sub>0</sub></td><td>239.252</td><td>82.1</td></tr><tr><td>H<sub>2</sub>CS 6<sub>4,3</sub>–5<sub>4,2</sub></td><td>206.002</td><td>245.1</td><td>CH<sub>3</sub>CCH 14<sub>1</sub>–13<sub>1</sub></td><td>239.248</td><td>93.3</td></tr><tr><td>H<sub>2</sub>CS 6<sub>4,2</sub>–5<sub>4,1</sub></td><td>206.002</td><td>245.1</td><td>CH<sub>3</sub>CCH 14<sub>2</sub>–13<sub>2</sub></td><td>239.234</td><td>114.7</td></tr><tr><td>SO 4<sub>5</sub>–3<sub>4</sub></td><td>206.176</td><td>38.6</td><td>CH<sub>3</sub>CCH 14<sub>3</sub>–13<sub>3</sub></td><td>239.211</td><td>150.3</td></tr><tr><td>SO<sub>2</sub> 3<sub>2,2</sub>–2<sub>1,1</sub></td><td>208.700</td><td>15.3</td><td>CH<sub>3</sub>CCH 14<sub>4</sub>–13<sub>4</sub></td><td>239.179</td><td>200.3</td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">(RxB lower sideband)</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>0,5</sub>–4<sub>0,4</sub><i>A</i></td><td>241.791</td><td>34.8</td><td>H<sub>2</sub>CS 7<sub>0,7</sub>–6<sub>0,6</sub></td><td>240.267</td><td>46.1</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>1,5</sub>–4<sub>1,4</sub><i>E</i></td><td>241.767</td><td>40.4</td><td>H<sub>2</sub>CS 7<sub>2,6</sub>–6<sub>2,5</sub></td><td>240.382</td><td>98.8</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>0,5</sub>–4<sub>0,4</sub><i>E</i></td><td>241.700</td><td>47.9</td><td>H<sub>2</sub>CS 7<sub>3,5</sub>–6<sub>3,4</sub></td><td>240.392</td><td>164.6</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>-1,4</sub>–4<sub>-1,3</sub><i>E</i></td><td>241.879</td><td>55.9</td><td>H<sub>2</sub>CS 7<sub>3,4</sub>–6<sub>3,3</sub></td><td>240.332</td><td>164.6</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>-2,4</sub>–4<sub>-2,3</sub><i>E</i></td><td>241.905</td><td>57.1</td><td>H<sub>2</sub>CS 7<sub>4,4</sub>–6<sub>4,3</sub></td><td>240.332</td><td>256.5</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>2,3</sub>–4<sub>2,3</sub><i>E</i></td><td>241.904</td><td>60.7</td><td>SO<sub>2</sub> 14<sub>0,14</sub>–13<sub>1,13</sub></td><td>244.254</td><td>93.9</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>2,4</sub>–4<sub>2,3</sub><i>A</i></td><td>241.842</td><td>72.5</td><td colspan="3" style="text-align:center; font-weight:bold;">(RxB upper sideband)</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>2,3</sub>–4<sub>2,2</sub><i>A</i></td><td>241.887</td><td>72.5</td><td>C<sub>2</sub>H <i>N</i> = 3–2, <i>J</i> = 7/2–5/2, <i>F</i> = 3–2</td><td>262.01</td><td>25.1</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>-3,3</sub>–4<sub>-3,2</sub><i>E</i></td><td>241.844</td><td>82.5</td><td>C<sub>2</sub>H <i>N</i> = 3–2, <i>J</i> = 7/2–5/2, <i>F</i> = 4–3</td><td>262.00</td><td>25.1</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>3,3</sub>–4<sub>3,2</sub><i>A</i></td><td>241.833</td><td>84.7</td><td>C<sub>2</sub>H <i>N</i> = 3–2, <i>J</i> = 5/2–3/2, <i>F</i> = 3–2</td><td>262.06</td><td>25.1</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>3,2</sub>–4<sub>3,1</sub><i>A</i></td><td>241.833</td><td>84.7</td><td>C<sub>2</sub>H <i>N</i> = 3–2, <i>J</i> = 5/2–3/2, <i>F</i> = 2–2</td><td>262.08</td><td>25.2</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>3,3</sub>–4<sub>3,2</sub><i>E</i></td><td>241.852</td><td>97.5</td><td>C<sub>2</sub>H <i>N</i> = 3–2, <i>J</i> = 5/2–3/2, <i>F</i> = 2–1</td><td>262.07</td><td>25.2</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>3,3</sub>–4<sub>3,2</sub><i>A</i></td><td>241.807</td><td>115.2</td><td>CH<sub>3</sub>CCH 15<sub>0</sub>–14<sub>0</sub></td><td>256.337</td><td>98.4</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>4,1</sub>–4<sub>4,0</sub><i>A</i></td><td>241.807</td><td>115.2</td><td>CH<sub>3</sub>CCH 15<sub>1</sub>–14<sub>1</sub></td><td>256.331</td><td>105.6</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>4,2</sub>–4<sub>4,1</sub><i>E</i></td><td>241.813</td><td>122.7</td><td>CH<sub>3</sub>CCH 15<sub>2</sub>–14<sub>2</sub></td><td>256.317</td><td>127.0</td></tr><tr><td>CH<sub>3</sub>OH 5<sub>-4,1</sub>–4<sub>-4,0</sub><i>E</i></td><td>241.830</td><td>130.8</td><td>CH<sub>3</sub>CCH 15<sub>3</sub>–14<sub>3</sub></td><td>256.292</td><td>162.7</td></tr><tr><td></td><td></td><td></td><td>CH<sub>3</sub>CCH 15<sub>4</sub>–14<sub>4</sub></td><td>256.258</td><td>212.6</td></tr><tr><td></td><td></td><td></td><td>H<sup>13</sup>CO<sup>+</sup> 3–2</td><td>260.255</td><td>25.0</td></tr><tr><td></td><td></td><td></td><td>SO 6<sub>6</sub>–5<sub>5</sub></td><td>258.256</td><td>56.5</td></tr><tr><td></td><td></td><td></td><td>SO 7<sub>6</sub>–6<sub>5</sub></td><td>261.843</td><td>47.6</td></tr></tbody></table>

Notes. The information is from the CDMS database (Müller et al. 2001).

clumps having constant accretion rates of $10^{-5} M_{\odot}yr^{-1}$ (dotted), $10^{-4} M_{\odot}yr^{-1}$ (dashed) and $10^{-3} M_{\odot}yr^{-1}$ (solid), for the most massive star in the cluster are also indicated (green lines). In these tracks, the other stellar members follow an equal accretion stopping probability, with an accretion rate $\propto M^{1.5}$. The orange pluses mark the time epoch of $2 \times 10^4$ yr for each accretion track. The SFE is assumed to be 30% and the underlying stellar population follows canonical IMF (Kroupa et al. 1993). Accretion luminosities are estimated by interpolating massive protostar models in Hosokawa & Omukai (2009).

The lower mass limit of $300 M_{\odot}$ corresponds to the mass of a massive clump in which at least one $>8 M_{\odot}$ star will form according to the normal IMF with an assumed star formation efficiency (SFE) of 30% (Kroupa et al. 1993; Sanhueza et al. 2017). The peak flux density of 2 Jy beam$^{-1}$ with beam full width at half maximum (FWHM) ~20" at 870 µm from ATLASGAL,

considering a distance of 6 kpc, implies a mass of $>100 M_{\odot}$ concentrated in the clump central ~0.5 pc region, assuming a dust temperature of 50 K and dust opacity of $1.8 \text{ cm}^2 \text{ g}^{-1}$ with a gas-to-dust ratio of 100. Thus, this selection criterion yields a sample of eight massive clumps (Table 2) with moderate to high central concentration. We note that our target sample is representative of more evolved sources with respect to those that fulfill these criteria. These sources can be easily detected with 1 mm lines given their favorable excitation conditions, and we further complement the sample with an infrared dark source, G18.

## 2.2. SMA observations

We performed SMA observations in the ~1 mm band toward seven clumps in the subcompact array configuration in June
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schematic picture of molecular cloud structure over spatial scales from >10 pc to ~0.1 pc: from cloud to core scale. A molecular cloud starts contraction from an initial stage that appears to be an infrared dark cloud (IRDC) and evolves into a star-forming one, embedding a number of molecular clumps. The massive clumps, ~1 pc in size, are generally composed of filamentary structures and cores at different evolutionary stages. In all figures, the yellow curved arrows represent turbulent motions and the purple arrows indicate gravitational contraction or gas inflows (along filaments). In the rightmost figure, the thick lines show filaments and the blue ovals indicate cores of different masses; the color gradient of the clump indicates a density gradient of the bulk gas. The characteristics of different structures are linked to textboxes by dotted arrows.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Luminosity–mass diagram of target sources (three-branched triangles). The evolutionary tracks of massive clumps of different envelope masses that form a cluster of stars with different accretion rates are shown as dotted and solid green lines; the evolutionary tracks of clumps of different envelope masses that are assumed to form a single massive star are shown as dash-dotted gray lines with arrows (Molinari et al. 2008; for more details see Sect. 2.1).</figcaption></figure>

2017, and in the compact array configuration in August 2017. The selected target sources are summarized in Table 2. Detailed information about each target source from previous studies is summarized in Appendix A. We used the dual receivers mode supported by the SMA Wideband Astronomical ROACH2 Machine (SWARM) backend: The RxA receivers covered the

frequency ranges of 188.4–196.7 GHz and 204.4–212.7 GHz in the lower and upper sidebands, respectively; the RxB receivers covered the frequency ranges of 238.5–246.8 GHz and 254.5–262.8 GHz, respectively. The intrinsic spectral channel width was 140 kHz. The molecular line transitions we covered are summarized in Table 1.

In addition, we retrieved archival SMA data toward the luminous ultra-compact HII region G10.624–0.38, which remains deeply embedded in a $M_{\text{gas}} = 10^3–10^4 M_{\odot}$ molecular clump and harbors a cluster of OB stars. These observations covered the CH$_3$OH $J = 5–4$ and $J = 7–6$ and the CH$_3$CN $J = 19–18$ line multiplets. We refer to Liu et al. (2010, 2011, 2012) for details of these observations.

We followed the standard SMA data calibration strategy. The application of system temperature ($T_{sys}$) information and the absolute flux, passband, and gain calibrations were carried out using the MIRIDL software package (Qi 2003). The absolute flux scalings were derived by comparing the visibility amplitudes of the gain calibrators with those of the absolute flux standard sources of the SMA, which were Callisto and Neptune for the subcompact and compact array observations, respectively. After calibration, we performed zeroth-order fitting of continuum levels from line-free channels and the joint weighted imaging adopting robust weighting of all continuum data were performed using the MIRIAD software package (Sault et al. 1995). The resultant synthesized beam is typically 4.′′5 at 241 GHz. The sensitivity ($3\sigma$) of continuum observation is ~0.04 Jy beam$^{-1}$ and of lines ~0.5 K. For clump G18 we do not obtain robust detection of the thermometer lines of CH$_3$CCH, H$_2$CS, and CH$_3$CN (Table 1) with this achieved line sensitivity. We used the previously published result from IRAM 30m telescope observations of 3 mm CH$_3$CCH and CH$_3$CN lines (Giannetti et al. 2017) instead.
<table><caption>**Table 2.** Target sources.</caption><thead><tr><th>Source <sup>(a)</sup></th><th>RA<br>(J2000)</th><th>Dec<br>(J2000)</th><th>Distance <sup>(b)</sup><br>(kpc)</th><th>Gas mass <sup>(c)</sup><br>(10<sup>2</sup> M<sub>☉</sub>)</th><th>Luminosity <sup>(c)</sup><br>(10<sup>3</sup> L<sub>☉</sub>)</th><th>L/M<br>(L<sub>☉</sub>/M<sub>☉</sub>)</th><th>Category</th></tr></thead><tbody><tr><td>G18.606–00.074 (G18)</td><td>18<sup>h</sup>25<sup>s</sup>08<sup>s</sup>.27</td><td>−12°45′22′′.7</td><td>3.7</td><td>19</td><td>0.75</td><td>0.4</td><td>IR weak</td></tr><tr><td>G28.397+00.080 (G28)</td><td>18<sup>h</sup>42<sup>s</sup>52<sup>s</sup>.08</td><td>−03°59′53′′.7</td><td>4.8</td><td>37</td><td>12</td><td>3.2</td><td>IR bright</td></tr><tr><td>G19.882–00.534 (G19)</td><td>18<sup>h</sup>29<sup>s</sup>14<sup>s</sup>.19</td><td>−11°50′28′′.4</td><td>3.7</td><td>20</td><td>99</td><td>4.8</td><td>IR bright</td></tr><tr><td>G08.684–00.367 (G08a)</td><td>18<sup>h</sup>06<sup>s</sup>23<sup>s</sup>.35</td><td>−21°37′05′′.2</td><td>4.8</td><td>27</td><td>27</td><td>10.2</td><td>IR weak</td></tr><tr><td>G31.412+00.307 (G31)</td><td>18<sup>h</sup>47<sup>s</sup>34<sup>s</sup>.32</td><td>−01°12′45′′.5</td><td>7.9</td><td>182</td><td>161</td><td>9.0</td><td>HII</td></tr><tr><td>G08.671–00.356 (G08b)</td><td>18<sup>h</sup>06<sup>s</sup>19<sup>s</sup>.23</td><td>−21°37′26′′.8</td><td>4.8</td><td>52</td><td>11</td><td>20.7</td><td>HII</td></tr><tr><td>G13.658–00.599 (G13)</td><td>18<sup>h</sup>17<sup>s</sup>23<sup>s</sup>.46</td><td>−17°22′09′′.2</td><td>4.5</td><td>12</td><td>24</td><td>19.9</td><td>IR bright</td></tr><tr><td>G10.624–00.380 (G10)</td><td>18<sup>h</sup>10<sup>s</sup>28<sup>s</sup>.638</td><td>−19°55′49′′.1</td><td>4.95</td><td>84</td><td>511</td><td>60</td><td>HII</td></tr></tbody></table>

Notes. <sup>(a)</sup>The abbreviated source names are given in parentheses. <sup>(b)</sup>The references for the source distances are: Green & McClure-Griffiths (2011) for G08a, G08b, and G19; Sanna et al. (2014) for G10; Wienen et al. (2015) for G13; Churchwell et al. (1990) for G31; Carey et al. (2000) for G28; Sridharan et al. (2005) for G18. <sup>(c)</sup>The gas mass and luminosity of the clump is taken from Urquhart et al. (2018), but scaled to the adopted distance of the clump.

## 2.3. APEX observations

Single-dish observations at 1 mm toward our target sources were performed with the MPIfR principal investigator (PI) instrument PI230 on the Atacama Pathfinder Experiment 12-m submillimeter telescope (APEX, Güsten et al. 2006), between March and September 2017 and July 2018 (Project M-099.F-9513A-2017, PI: Yuxin Lin). The PI230 receiver is a dual polarization sideband separating heterodyne system with a total of 32 GHz bandwidth working at 230 GHz, and can cover the spectral range of 200–270 GHz. We conducted observations with two spectral setups, covering frequency ranges of 202.2–210.0 GHz, 218.0–225.8 GHz and 239.2–247.0 GHz, 255.0–262.8 GHz. For each target source, a region of 3' × 3' centered at the source position was mapped in the on-the-fly (OTF) mode with both setups.

During the observations, the typical precipitable water vapor (PWV) was ~1.5–2.5 mm. The pointing was determined by continuum observations on Saturn when available, or CO $J = 2-1$ observations on bright nearby evolved stars (e.g., RAFGL2135, R-Dor). The pointing accuracy was found to be within 3''. The focus was checked on Saturn every 2–4 h. The main beam efficiency ($\eta_{mb}$) for the PI230 instrument varied over the observing period, with a range of ~63–72%<fnref n="1" />. The calibration uncertainty was typically within 20%, estimated based on the flux measurement of the pointing sources.

Basic data reductions were done with the GILDAS software package<fnref n="2" />, including flagging of bad spectra, baseline subtraction, unit conversion ($T_A^*$ to $T_{mb}$), and building spectral cubes. Final spectral cubes were re-sampled to 0.5 km s⁻<fnref n="1" /> spectral resolution for all lines.

## 2.4. SMA-APEX combination and imaging

For our primary target lines covered by both SMA and APEX, which have extended emission, namely CH₃CCH $J = 12-11$, CH₃OH $J = 5-4$, C₂H $J = 3-2$, H¹³CO⁺ $J = 3-2$, CS $J = 5-4$, C³⁴S $J = 5-4$, and SO $J = 4_5-3_4$, we combined the two data sets in the Fourier domain (uv-domain) with MIRIAD. This combination is essentially imaging together the pseudo-visibilities

generated from single-dish observations and interferometer measurements, so that the short-spacing information including zero baseline that is obtained with a single dish can be complemented by the interferometry data; the method is commonly referred to as the joint deconvolution or joint reconstruction method (Kurono et al. 2009; Koda et al. 2011).

In the combination procedure, we first deconvolved the APEX data from the single-dish Gaussian beam (FWHM ~30'') and then multiplied the resultant image with the primary beam of the SMA observations. The obtained image was then used to generate a single-dish uv model (i.e., the pseudo-visibilities) by randomly sampling a visibility distribution to match that of the single-dish beam pattern. The zero spacing visibility was also sampled and added to the produced pseudo-visibilities. Finally, the pseudo- and interferometric visibilities were imaged together to produce the combined image. In the final imaging step, we applied a Gaussian tapering function of FWHM ~2'' to increase the detectability of extended emission. In the end we adopted a final step to linearly combine the product of this standard joint deconvolution method with the single-dish image in the Fourier domain, using immerge in the MIRIAD package. This step was necessary and was found to preserve the single-dish overall fluxes better than using solely the joint deconvolution method, due to the fact that the deconvolution method is not flux conserving; a similar procedure is adopted in Monsch et al. (2018), among others. The combined images have comparable total fluxes to the APEX data, within a difference of 15%.

For details on the proper weighting scheme in the joint deconvolution method and the impact of sensitivities of single-dish and interferometry data, we refer to Kurono et al. (2009).

## 2.5. Ancillary data: mid- to far-infrared and submm continuum data from multiple single-dish telescopes

We used the single-dish mid- and far-infrared, and submillimeter (submm) continuum data to constrain the bulk gas structures and construct the SEDs (Fig. 3). In addition to the 870 µm data from ATLASGAL survey (Schuller et al. 2009; Csengeri et al. 2016) obtained by APEX-LABOCA (Siringo et al. 2009), we also adopted 350 µm data obtained by the CSO-SHARC2 or APEX-SABOCA instrument. The information of the observations and data reduction procedure are detailed in Lin et al. (2017) and Lin et al. (2019). For sources without available 350 µm from ground-based telescope (of 10'' angular resolution), we used

<footnote n="1">¹ http://www.apex-telescope.org/telescope/efficiency/</footnote>

<footnote n="2">² https://www.iram.fr/IRAMFR/GILDAS/</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Overall work flow showing the radiative transfer modeling procedure followed in this work.</figcaption></figure>

the available observations from *James Clerk Maxwell Telescope* (JCMT)<fnref n="3" /> Submillimetre Common-User Bolometer Array 2 (SCUBA2; Dempsey et al. 2013; Chapin et al. 2013) at 450 µm from the online data archive.

We also retrieved the archival *Herschel*<fnref n="4" /> images at 70/160 µm and 250/350/500 µm from the *Herschel* Infrared Galactic Plane (Hi-GAL) survey (Molinari et al. 2010) taken by the PACS (Poglitsch et al. 2010) and SPIRE instruments (Griffin et al. 2010). For the mid-infrared data, we used the 24 µm images

from the MIPS Galactic Plane Survey (MIPSGAL, Carey et al. 2009) taken by *Spitzer* telescope.

# 3. Results and analysis

## 3.1. *Outline of the modeling and analysis procedure*

In this work, we aim to provide a description of the gas density and temperature of massive clumps by performing radiative transfer calculations of molecular lines and the multi-wavelength dust continuum. In this section, we introduce the workflow of the modeling, starting by basic definitions of molecular line excitation. The modeling steps (shown in Fig. 3) are explained in more detail in Sects. 3.4–3.6 and in Appendices E and F. The results of the radiative transfer models are discussed in Sect. 4.

Massive clumps have average molecular hydrogen densities of typically ~$10^4$ cm⁻<fnref n="3" /> (Csengeri et al. 2014; Urquhart et al. 2018); the collisional partner participating in the excitation and de-excitation of molecular lines considered in this paper is primarily hydrogen gas. The critical density ($n_{crit}$) (Table 3) defines the way in which a molecule in an excited state decays to ground state. When the main collisional partner is hydrogen, it stands for the critical hydrogen density at which timescales of radiative decay and collisional de-excitation are comparable; in other words, the net radiative decay rate from level *J* to a certain lower level equals the rate of collisional de-population out of the upper level *J* for a multi-level system (e.g., Wilson et al. 2013).

<footnote n="3">³ The *James Clerk Maxwell* Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The *James Clerk Maxwell* Telescope has historically been operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada and the Netherlands Organization for Scientific Research. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation.</footnote>

<footnote n="4">⁴ *Herschel* is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>
<table><caption>**Table 3.** Critical density for transitions of interest.</caption><thead><tr><th>Transitions</th><th>Critical density <sup>(a)</sup><br>(cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>CS 5-4</td><td>1.3 × 10<sup>6</sup></td></tr><tr><td>CH<sub>3</sub>CCH 12<sub>0</sub>–11<sub>0</sub></td><td>2.9 × 10<sup>4</sup> <sup>(b)</sup></td></tr><tr><td>H<sub>2</sub>CS 6<sub>0.6</sub>–5<sub>0.5</sub></td><td>2.2 × 10<sup>5</sup></td></tr><tr><td>CH<sub>3</sub>OH 5<sub>1.5</sub>–4<sub>1.4</sub></td><td>6.3 × 10<sup>5</sup></td></tr><tr><td>CH<sub>3</sub>OH 5<sub>4.5</sub>–4<sub>4.4</sub></td><td>5.3 × 10<sup>7</sup></td></tr><tr><td>CH<sub>3</sub>CN 13<sub>0</sub>–12<sub>0</sub></td><td>4.6 × 10<sup>6</sup></td></tr><tr><td>CH<sub>3</sub>CCH 14<sub>0</sub>–13<sub>0</sub></td><td>4.3 × 10<sup>4</sup> <sup>(b)</sup></td></tr><tr><td>H<sub>2</sub>CS 7<sub>0.7</sub>–6<sub>0.6</sub></td><td>3.7 × 10<sup>5</sup></td></tr><tr><td>C<sub>2</sub>H N = 3–2, J = 7/2–5/2, F = 4–3</td><td>2.7 × 10<sup>5</sup></td></tr><tr><td>SO 4<sub>5</sub>–3<sub>4</sub></td><td>2.2 × 10<sup>5</sup></td></tr><tr><td>SO 6<sub>6</sub>–5<sub>5</sub></td><td>7.0 × 10<sup>5</sup></td></tr><tr><td>SO 7<sub>6</sub>–6<sub>5</sub></td><td>1.0 × 10<sup>6</sup></td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup> 3–2</td><td>9.5 × 10<sup>5</sup></td></tr><tr><td>SO<sub>2</sub> 14<sub>0.14</sub>–13<sub>1.13</sub></td><td>4.0 × 10<sup>5</sup></td></tr><tr><td>SO<sub>2</sub> 3<sub>2.2</sub>–2<sub>1.1</sub></td><td>1.7 × 10<sup>5</sup></td></tr><tr><td>CH<sub>3</sub>CCH 15<sub>0</sub>–14<sub>0</sub></td><td>5.7 × 10<sup>4</sup> <sup>(b)</sup></td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Calculated following the definition in Shirley (2015) in the optically thin limit ($\tau \sim 0.1$) at 50 K, considering a multi-level energy system whenever necessary. <sup>(b)</sup>Calculated using collisional coefficients of CH<sub>3</sub>CN.

With gas densities close to and well above $n_{crit}$, the thermalization of energy levels is achieved, such that the excitation temperature ($T_{ex}$) can approximate the gas kinetic temperature ($T_{kin}$), with the population of energy levels reaching Boltzmann prediction (local thermodynamic equilibrium, LTE). On the other hand, if gas densities are below $n_{crit}$ (subthermal excitation), then the population of the upper energy level is sensitive to varying gas densities. Observations of multiple transitions with different $n_{crit}$ can probe a range of gas densities by showing rather different ratios of line intensities. In particular, if the energy levels associated with these transitions are of similar energy, then the dependence of line ratios on temperature is minimized, and so does the degeneracy of the mutual effect of gas temperature and density in determining level populations. We take advantage of these radiative properties to use selected molecular transitions as densitometers of our target sources.

In a simple view, massive star-forming clumps may be thought of as multi-layered gas structures showing centrally peaked gas density profiles. This is a natural outcome under self-gravity. From the outermost layer to the innermost region, transitions of higher and higher $n_{crit}$ are thermalized progressively. Using a combination of thermometers of different $n_{crit}$, based on LTE assumption, can constrain gas temperatures over continuous spatial scales (with respect to the clump center). Analogous to a densitometer, a thermometer is defined here as a set of molecular lines of a certain species whose level population is only (or dominantly) sensitive to gas temperature, which arise from energy levels spanning a wide energy range and are connected ideally only through collisions, provided for example by K-ladder lines of symmetric top molecules.

Considering the gas density regime of massive star-forming clumps, and based on previous single-dish experiments (Giannetti et al. 2017), we have identified CH<sub>3</sub>CCH, H<sub>2</sub>CS, and

CH<sub>3</sub>CN lines at the 1 mm band (listed in Table 1) as ideal tracers for measuring temperature profiles of massive clumps. These tracers have $n_{crit}$ of several $10^4$, $10^5$, and $10^6$ cm<sup>-3</sup>, respectively (Table 3). On the other hand, the combination of distinct $n_{crit}$ triggers a filtering effect, such that with each thermometer the region it probes is limited to the gas density regimes ranging around and above its $n_{crit}$. Contamination by foreground and background gas layers of lower density to the observed emission is therefore negligible. This means that the line-of-sight (averaging) effect is reduced to the gas component of a limited density range.

In addition, optical depths ($\tau$) are low when typical abundances and excitation conditions are considered for these molecules, which implies that line ratios probe the gas kinetic temperature at the inner location of the gas layer. With these properties in mind, we derive the rotational temperature ($T_{rot}$) maps under the LTE assumption using multiple thermometers in Sect. 3.4. Temperature measurements of the outer regions are obtained using the extended CH<sub>3</sub>CCH and H<sub>2</sub>CS emission, and combined with temperatures derived from CH<sub>3</sub>CN, which is confined to the central region of the clumps. This combination allows us to establish the full radial temperature profile of the clumps. We also use multi-wavelength single-dish dust emission (SD continuum, as in Fig. 3) to derive dust temperature maps by building spectral energy distributions (SEDs) assuming single-component modified blackbody emission (Lin et al. 2016, 2019). The dust temperatures at the outer layer of clumps, are used to complete the temperature profile at larger radii for the clumps. With the simple one-component LTE modeling and one-component dust SED construction, we derive the projected radial profile of the obtained multiple temperature maps as the radial profiles, denoted $T(r)$. With this approximation, a natural difference caused by line-of-sight (LOS) projection effects may appear as a function of radius due to density-weighted emission. However, as previously mentioned, due to the density-filtering effect by combination of multiple thermometers, the difference between the two profiles is largely minimized. Moreover, the projected radial temperature profile used as radial temperature profile is further benchmarked and refined by SED construction from a full radiative transfer modeling of dust based on a density profile adopted for the clump (Sect. 3.6, Appendix E), and further shown to be able to produce the observed CH<sub>3</sub>CCH lines and their spatial variation by full line radiative transfer models (Appendix F).

To probe the gas density, we rely on the CH<sub>3</sub>OH line series in the 1 mm band as a densitometer. We adopt one-component non-LTE models to derive the hydrogen volume density ($n(H_2)$) maps and benchmark the results using full non-LTE radiative transfer modeling (Sect. 3.6, Appendix F). The highest and lowest $n_{crit}$ of the adopted line series are listed in Table 3. Moreover, with measured radial temperature profiles, the degeneracy of temperature and density can be further reduced by introducing $T_{kin}$ ($T(r)$) in the non-LTE modeling, to constrain solely $n(H_2)$. We adopt this strategy in Sect. 3.5 (see also Appendix C). In Sect. 3.5, we use the one-component non-LTE model of CH<sub>3</sub>OH lines to constrain $n(H_2)$. As mentioned before, the excitation of molecular lines has a selective effect on gas densities. Since we are also interested in the bulk gas density structure of massive clumps, we complement the density distribution measure with single-dish multi-wavelength dust emission. We conduct full radiative transfer modeling of dust to fit these data (Appendix E), again incorporating the $T(r)$ initially measured from thermometers. To benchmark the one-component non-LTE models, and to
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spitzer IRAC RGBs (R: 8.0 µm; G: 4.5 µm; B: 3.6 µm) maps of the target sources. Yellow contours show the ATLASGAL 870 µm emission from 1 Jy beam⁻¹ to the peak flux for each source, using seven levels with uniform spacing. White contours show SMA 1.2 mm emission from 3σ to the peak flux using five levels with uniform spacing. Negative flux levels of the 1.2 mm continuum are shown in contours of dotted lines, from −1σ to the minimum negative flux with two levels. The beam of the SMA continuum is shown in the lower left corner of each plot. The beam size of the archival data for source G10 is much smaller than other sources (Sect. 2.2). The primary beam size is shown in each plot as a blue dashed circle.</figcaption></figure>

understand the difference of gas density results between modeling dust emission and simple non-LTE modeling of CH₃OH, we utilize non-LTE full radiative transfer calculation of lines in Appendix F, for CH₃OH and CH₃CCH lines. This also helps us to examine the possibility of spatial abundance variations of these lines as an additional factor in affecting the distribution of line emission. In this effort, particularly, the full non-LTE modeling of CH₃CCH provides a sanity check on the measured radial temperature profile *T(r)* from rotational temperature maps. The workflow of the whole procedure is graphically summarized in Fig. 3.

## 3.2. SMA 1.2 mm continuum

The SMA 1.2 mm dust continuum images resolved two compact sources (separated by ~7.''2, ~0.15 pc) in G18, and resolved isolated compact sources in the rest of the samples (Fig. 4). Hereafter we refer to these compact sources as core structures. Before any further analyses, we utilized the archival centimeter band data to subtract free-free contamination in G08b, G31, and G10, assuming optically thin emission (i.e., $S_v^{\text{free-free}} \propto v^{-0.1}$, for details see Appendix B). Then the core radius is defined as the area above 5σ emission contours of the 1.2 mm images. The core effective radius, peak flux, and integrated flux are listed in Table 4.

We assumed that dust emission in all cores is optically thin at 1.2 mm. Based on the OH5 opacity model (i.e., $\kappa_{1.2\text{mm}} = 0.81 \text{ cm}^2 \text{ g}^{-1}$; Ossenkopf & Henning 1994), we converted the continuum intensity detected at >5σ to dust mass surface density, which was subsequently converted to gas mass surface density by assuming that the gas-to-dust mass ratio is 100. In these mass estimates, we assumed that dust temperature is identical to the gas temperature *T(r)* which we derived (and refined, Eq. (3)) (see Sects. 3.4, 3.6, and Appendix E).

<table><caption>**Table 4.** Source properties from 1.2 mm SMA continuum.</caption><thead><tr><th>Source</th><th>S<sub>1.2mm</sub><br>(Jy)<sup>(a)</sup></th><th>R<sub>eff</sub><br>(pc)<sup>(b)</sup></th><th>T&#773;<sub>d</sub><br>(K)<sup>(c)</sup></th><th>M<sub>core</sub><br>(M<sub>☉</sub>)</th><th>M<sub>core</sub><sup>Abel</sup><br>(M<sub>☉</sub>)</th><th>ρ&#773;<sub>core</sub><br>(10<sup>4</sup> cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>G18<sup>(e)</sup></td><td>0.13<br>0.02</td><td>0.10<br>0.04</td><td>36<br>19</td><td>24.0<br>5.5</td><td>22.1<br>—</td><td>6.0<br>—</td></tr><tr><td>G28</td><td>0.69</td><td>0.12</td><td>49</td><td>93.5</td><td>103.7</td><td>19.5</td></tr><tr><td>G19</td><td>1.07</td><td>0.10</td><td>37</td><td>124.1</td><td>65.3</td><td>23.6</td></tr><tr><td>G08a</td><td>0.35</td><td>0.10</td><td>50</td><td>49.1</td><td>18.6</td><td>8.0</td></tr><tr><td>G31</td><td>5.12</td><td>0.17</td><td>94</td><td>943.8</td><td>914.0</td><td>61.0</td></tr><tr><td>G08b</td><td>1.91</td><td>0.16</td><td>56</td><td>223.2</td><td>188.1</td><td>15.0</td></tr><tr><td>G13</td><td>0.51</td><td>0.11</td><td>48</td><td>59.0</td><td>58.5</td><td>15.2</td></tr><tr><td>G10</td><td>1.85</td><td>0.04</td><td>132</td><td>102.0</td><td>67.8</td><td>210.0</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Total flux above 5σ. For source G08b, G31 and G10, subtraction of free-free emission is considered (details in Appendix B). <sup>(b)</sup>Effective radius is defined as $\pi R_{\text{eff}}^2 = \text{Area}$, where Area stands for the emission region above 5σ for each 1.2 mm map, i.e., non-deconvolved averaged size. <sup>(c)</sup>Dust temperature is assumed to be equal to gas temperature *T(r)* obtained and refined in Sect. 3.4, and the average temperature $\bar{T}_d$ is calculated by averaging over pixels that have continuum emission above 5σ. <sup>(d)</sup>Average density calculated from $M_{\text{core}}^{\text{Abel}}$ and $R_{\text{eff}}$. <sup>(e)</sup>For clump G18, the calculation for two cores are given. In the Abel inversion calculation which assumes spherical symmetry, only the central, more massive core is considered. So $M_{\text{core}}^{\text{Abel}}$ and $\bar{\rho}_{\text{core}}$ are omitted for the secondary core.

There is a subtlety in the way we applied *T(r)*, which is related to the assumption of the thermal and density structures of the cores. We compared two ways of applying *T(r)*. In the first, we defined a mean core gas temperature $T_{\text{core}}$ by making averages of *T(r)* within the core size. For each pixel, we then adopted a dust temperature which is equal to $\min\{T(r = \ell), T_{\text{core}}\}$ when deriving $n(H_2)$, where $\ell$ is the projected distance from the
pixel to the 1.2 mm continuum peak (i.e., centers of the sources). This means that we use the lower value of the two temperatures (average core temperature and radial temperature) at each pixel position to estimate the gas mass probed by dust emission. This is a reasonable assumption since dust emission is sampling all the gas components along the LOS and the average mass temperature is likely dominated by the outer, colder gas component. Given that the projected radius $l$ is always smaller than the radius $r$, this approach still tends to overestimate the dust temperatures at small projected radii, although it is alleviated. This in turn results in an underestimate of $n(H_2)$.

In a second approach, we assumed that the cores are spherically symmetric and optically thin. We used an Abel transform to convert the observed azimuthally averaged intensity profile of 1.2 mm emission to gas density $\rho(r)$ (for more details, see Roy et al. 2014), as

$$ \rho(r) = -\frac{1}{\pi \kappa_v B_v [T(r)]} \int_r^{r_{\text{eff}}} \frac{dI_v}{db} \frac{db}{\sqrt{b^2 - r^2}}, \quad (1) $$

where $r_{\text{eff}}$ is the core effective radius. We then integrated $\rho(r)$ over the LOS to obtain another version of the $n(H_2)$ map. The two versions of the $n(H_2)$ maps agree within a factor of 1.5–2. The average dust temperature within the core, two sets of mass estimates $M_{\text{core}}$ and $M_{\text{Abel}}$, and average core density are summarized in Table 4.

### 3.3. Distribution of emission from CH$_3$CCH, H$_2$CS, CH$_3$CN, CH$_3$OH lines and 1.2 mm continuum

Figure 4 shows the 1.2 mm dust continuum images taken with the SMA. We resolved two compact sources (separated by ~7.7′′2, ~0.15 pc) in G18, and resolved isolated compact sources in the rest of the samples. Figures 5–6 show the integrated intensity maps of CH$_3$CCH, H$_2$CS, and CH$_3$CN, which are overlaid on the integrated intensity maps of CH$_3$OH. In general, the CH$_3$OH lines and the lower K ladders of CH$_3$CCH were resolved on 0.3–0.4 pc scales, while the CH$_3$CN lines and higher K ladders of H$_2$CS were resolved on 0.1–0.2 pc scales. The results of our quantitative analyses are presented in the following subsections.

## 3.4. Deriving pixel-based gas rotational temperature maps with LTE modeling for multiple thermometers

### 3.4.1. Thermometers

CH$_3$CCH and CH$_3$CN are symmetric top molecules. Their K-ladder populations at a certain J level are determined primarily through collisions. Therefore, they have been regarded as thermometers for molecular clouds (Kuiper et al. 1984; Bergin et al. 1994). Given their similar geometry and molecular weight, CH$_3$CCH and CH$_3$CN are often assumed to have the same collisional coefficients, while CH$_3$CN has higher dipole moments than CH$_3$CCH. This means that a molecular clump can exhibit brighter CH$_3$CN line emission than CH$_3$CCH even when the excitation of the CH$_3$CN molecules is limited to small pockets of dense gas (e.g., confined to the hot core region).

The CH$_3$CN lines have been very commonly observed (Cummins et al. 1983; Sutton et al. 1986; Fayolle et al. 2015). They have been regarded as good tracers of hot molecular cores owing to the fact that they were mainly detected around significantly heated regions. On the other hand, CH$_3$CCH has been detected in spatially more extended, lower temperature regions (e.g., Bergin et al. 1994; Öberg et al. 2014) and is therefore

particularly advantageous for probing sources in relatively early evolutionary stages (Molinari et al. 2016), prior to hot core formation. Giannetti et al. (2017) showed that among the various thermometers the kinetic temperature constrained by CH$_3$CCH is representative of gas temperatures of massive clumps.

The H$_2$CS molecule is a near-prolate rotor. Its transitions between levels at various K ladders are also excellent indicators of the gas kinetic temperature (Blake et al. 1994). As a sulfur-bearing species, the gas phase H$_2$CS abundance can be enhanced either by direct evaporation or by outflow or shock sputtering (e.g., Bachiller & Gutiérrez 1997; Minh et al. 2011). Previous observations have also revealed that the H$_2$CS emission originates from extended warm regions surrounding compact hot cores (e.g., Helmich & van Dishoeck 1997).

### 3.4.2. Modeling procedure

Using the XCLASS package (Möller et al. 2017), we established a pixel-by-pixel LTE model-fitting procedure for the observed CH$_3$CCH, H$_2$CS, CH$_3$CN, and CH$_3$OH $v_t = 1$ lines (for G10), which returns the best fit of source size, rotational temperature ($T_{\text{rot}}$), molecular column density ($N_{\text{mol}}$), linewidth ($\Delta V$), and the source velocity ($V_{\text{source}}$). In this specific implementation, we fixed the source size to the synthesized beam size (i.e., assuming beam filling factor of 1) and optimized the rest of the free parameters. The optimization procedure employed an initial global parameter search using the bees algorithm (Pham et al. 2006), which was followed by Levenberg-Marquardt iterations. This procedure helps avoid trapping in local minima. We fit the $J = 15–14$ and $J = 14–13$ ladders of CH$_3$CCH together, and the $J = 12–11$ ladders separately, given that the former lines show less extended emission and appear to trace hotter gas. Examples of the fitted spectra are presented in Fig. 7. Examples of the obtained rotational temperature maps are shown in Fig. 8. When deriving these rotational temperature maps, we use pixels where the intensity of the third lowest energy transition in consideration is greater than our $3\sigma$ detection limit. For the hot molecular core G31, in addition to the few molecular species that we targeted, the fittings also considered several other species that can potentially make a prominent contribution in our spectra. They are shown in Fig. 7 for G31 in different colors. For clump G18 we do not have robust detection of these thermometer lines from our SMA observations, and we rely on previous IRAM 30 m telescope observations (Giannetti et al. 2017) to describe the temperature profile, with the same thermometers but their lower transitions. These pointed observations from the IRAM 30 m telescope did not give information on relevant physical scales; rather, based on a fixed temperature profile ($T(r) \propto R^{0.4}L^{-0.25}$, Giannetti et al. 2017) the radius of a certain measured rotational temperature was deduced. Following the same workflow (Fig. 3), these measurements are combined with one-component dust temperatures (Lin et al. 2019) in the outer region of the clump, to compose an initial temperature profile $T(r)$ to be refined later by SED comparison using RADMC-3D.

## 3.5. Deriving the pixel-based hydrogen volume density maps with non-LTE RADEX model

### 3.5.1. Methanol lines

Methanol (CH$_3$OH) is a slightly asymmetric top molecule. It has three types of symmetry, which are denoted A, E1, and E2. The E1 and E2 states can be thought of as doubly degenerate states of the E symmetry where the quantum number of the angular momentum along the symmetry axis of the CH$_3$ group ($k$) can
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Integrated intensity maps (gray contours) of CH₃CCH, H₂CS, and CH₃CN toward sources G19, G08a, and G08b. Integrated intensity of CH₃OH 5₀.₅–4₀.₄ ([$v_{lsr}$ − 3, $v_{lsr}$ + 3] km s⁻¹) is shown in color scale. Gray contours show the intensity levels with uniform intervals from 5σ up to the peak flux, with the emission range (Jy beam⁻¹ km s⁻¹) indicated in the lower left corner of each panel. Colored contour in the *left panel* shows the location of 0.8 × peak emission of the 1.2 mm SMA continuum image. The green ellipses indicate beams of corresponding molecular lines (open) and CH₃OH 5₀.₅–4₀.₄ line (filled).</figcaption></figure>

take either positive and negative values. The torsional ground state transitions E–CH₃OH 5<sub>k,5</sub> – 4<sub>k,4</sub>, K = 0, ±1, ±2, ±3, ±4 (v<sub>t</sub> = 0) were found to be a good densitometer for gas denser than 10<sup>4</sup> cm<sup>-3</sup> (Leurini et al. 2004, 2007).

The excitation of these *K* ladders is usually observed to be
subthermal. These *K* = 0 and *K* = ±1 ladders occupy a rather
narrow range of upper level energies (*E*<sub>up</sub> ~ 40–55 K). At the
same time, they cover a wide range of critical densities (~10<sup>5</sup> to
≥10<sup>7</sup> cm<sup>-3</sup>) (Table 3), which implies that the line ratios of two
*K* components can be good density probes. The higher *K* com-
ponents (*K* ≥ 3, *E*<sub>up</sub> > 80 K) are generally excited in hot regions

where the gas volume densities are close to or higher than the
critical densities. Hence, ratios of the *K* ≥ 3 components pro-
vide additional constraints on kinetic temperature. Apart from
the high abundance of CH₃OH, it is this property of the methanol
energy system and the relatively low upper level energies of the
*K* < 3 transitions that make this line series sensitive to gas den-
sity for a broad range of physical conditions in molecular clouds
(Leurini et al. 2004).

As illustrated in Fig. 5, the CH₃OH emission appears clumpy and exhibits elongated structures, extending for up to 0.5 pc with respect to the continuum peak. The *K* < 2 transitions of
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as Fig. 5, but for target clumps G13, G28, and G31.</figcaption></figure>

E–CH₃OH (5–4) (νᵢ = 0) are excited over an extended region, while the emission of the K > 2 lines are confined to the central region of the clumps.

### 3.5.2. Modeling procedure

We produced a series of large velocity gradient (LVG) RADEX models (van der Tak et al. 2007) to search for the best fits of *n*(*H*₂), CH₃OH column density, N(CH₃OH-*E*)/N(CH₃OH-*A*), and kinetic temperature (*T*<sub>kin</sub>) to the observed CH₃OH (*J* = 5–4, νᵢ = 0) lines. We took the collisional rates from Rabli & Flower (2010), which were evaluated for temperatures from 10 to 200 K. We adopted a Markov chain Monte Carlo (MCMC) method to derive the parameters and estimate the associated uncertainties, taking into consideration the upper limits for weakly detected line components. The details of the modeling procedure are elaborated in Appendix C where the formulas used for the likelihood function are given. In the fitting, for each pixel we enforced the

posterior distribution of $T_{\text{kin}}$ to be a narrow Gaussian distribution centralized at $T(r)$ (see Sect. 3.6), as measured in Sect. 3.4 from the multiple rotational temperature maps. Although the ratios between the lower *K* ladders of CH₃OH lines depend only weakly on the kinetic temperature, having a fixed term helps avoid randomly converged parameters, which is useful to ensure that the resultant parameter maps are continuous. The obtained *n*(*H*₂) maps are shown in Fig. 9. The CH₃OH column density maps are shown in Fig. C.2.

### 3.6. Radial density and temperature profiles used in full radiative transfer models

We used the RADMC-3D code (Dullemond et al. 2012) in our full radiative transfer analyses (Fig. 3; Sect. 3.1) for the multi-wavelength dust continuum (Appendix E). We assumed that the gas density profile for the bulk gas ($\rho_{\text{bulk}}(r)$) is described by the
<figure><img src="image_8.png" /><figcaption>Fig. 7. Example spectra of thermometer lines CH₃CN J = 13-12, CH₃CCH J = 14-13, H₂CS J = 6-5, H₂CS J = 7-6 at the continuum peak of the target source. The blue profiles show the XCLASS LTE fitting results. For source G31.412+0.307, which presents significant line blending from other species, the fittings also included those species and transitions that can potentially make prominent contributions to the spectrum.</figcaption></figure>

functional form

$$ \rho_{\text{bulk}}(r) = \begin{cases} \bar{\rho} \cdot \left(\frac{r}{r_c}\right)^q & (r \le R_{\text{clump}}), \\ 0 & (r > R_{\text{clump}}), \end{cases} \qquad (2) $$

where $\bar{\rho}$ is the mean hydrogen gas number density, and $r_c$ is the radius where $\rho(r) = \bar{\rho}$, $R_{\text{clump}}$ is the assumed outer radii of the clumps that were fixed to the FHWM measured from the ATLASGAL 870 µm maps (see Contreras et al. 2013). When converting gas density to mass density, we assume that the mass per hydrogen molecule is 2.8 $m_{\text{H}}$, where $m_{\text{H}}$ is the hydrogen atom mass. We assumed that the gas-to-dust mass ratio is 100.

We parameterized the measured temperature profiles $T(r)$ by

$$ T(r) = \omega T_{\text{in}} \left(\frac{r}{r_{\text{in}}}\right)^{-0.5} + (1 - \omega)T_{\text{out}}, \qquad (3) $$

where $\omega = e^{-r_{\text{out}}/r_c}$ is an exponential tapering function characterized by outer radius $r_{\text{out}}$, and $T_{\text{in}}$ and $T_{\text{out}}$ are the characteristic

temperatures at the radius $r_{\text{in}}$ and at asymptotically large radii, respectively. In this equation, the first term describes radiative heating by the centrally embedded stars, while the second term can be attributed to the ambient radiation fields of the massive clumps (Liu et al. 2019). The multiplicative factors $\omega$ and $(1 - \omega)$ prescribe the transition from one heating regime to the other.

Based on the multiple rotational temperature maps, we derived the (projected) radially averaged temperature profile and obtained best-fit parameters $T_{\text{in}}$, $T_{\text{out}}$, and $r_{\text{out}}$, while $r_{\text{in}}$ was initially kept as a fiducial value of 0.02 pc. Figures 10 and 11 show the comparison between these multi-wavelength radial intensity profiles and the SEDs evaluated from best-fit RADMC-3D (Dullemond et al. 2012) models. Based on this comparison, we re-adjusted the $r_{\text{in}}$ of Eq. (3) in the RADMC-3D modeling to obtain an SED profile consistent with the observed data points. The parameter $T(r)$ was then updated by the refined temperature profile. The parameters $T_{\text{in}}$, $T_{\text{out}}$, $r_{\text{in}}$, and $r_{\text{out}}$ that define $T(r)$ are
<figure><img src="image_9.png" /><figcaption>Fig. 7. continued.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Rotational temperature maps of source G08b derived from multiple thermometer lines using the XCLASS package (Sect. 3.4). The green contours indicate SMA 1.2 mm continuum levels from 0.3 to 0.9 × peak flux (Table 4) in five steps of uniform interval. The beams of the continuum and the respective lines are shown in the lower left corner as green and hatched ellipses, respectively.</figcaption></figure>

listed in Table 5. With $T(r)$ defined, we fixed the dust temperature profile in the multi-wavelength continuum modeling for the bulk gas and obtain $\rho_{bulk}(r)$. The best-fit parameters $\bar{\rho}$ and $q$ as in Eq. (2) are listed in Table 6.

From RADEX modeling of CH$_3$OH lines we constrained the radial density profiles for the dense gas, $\rho_{dense}(r)$, from the obtained $n(H_2)$ maps. Similarly, when deriving the $n(H_2)$ maps (Fig. 9) we fixed the gas kinetic temperature in the modeling to $T(r)$ for each pixel. We adopted a single power-law form as Eq. (2) to characterize the dense gas density profiles,

$$ \rho_{dense}(r) = \rho_{0.01\,pc} \left( \frac{r}{0.1\,pc} \right)^{q_{radex}}, \quad (4) $$

where $\rho_{0.01\,pc}$ is the reference gas density at 0.1 pc. The description is valid up to a maximum scale of $R_{max}$, which is determined from the largest radius where $n(H_2)$ can be robustly estimated. These parameters are also listed in Table 5. Figure 12 shows the comparison between the model fits and the observed radial profiles.

We then conducted full radiative transfer modeling with LIME (Brinch & Hogerheijde 2010) to benchmark and refine these results. In the LIME modeling of CH$_3$OH and CH$_3$CCH lines, we first adopted the bulk gas density profile $\rho_{bulk}(r)$ constrained from single-dish dust continuum modeling, and $T(r)$

with assumed abundance profiles to find the best-fit models. We parameterized the molecular abundance profiles ($X_{mol}(r)$) as

$$ X_{mol}(r) = \begin{cases} X_{\text{out}} & (T(r) < T_{\text{jump}}), \\ X_{\text{in}} \equiv f_{\text{inc}}X_{\text{out}} & (T(r) > T_{\text{jump}}), \end{cases} \quad (5) $$

where $T_{jump}$ is a threshold temperature chosen to be either 30 or 80 K, $X_{out}$ is the abundance at the outer radii, and $f_{inc}$ is an increment factor to characterize the abundance enhancement in the inner regions of higher temperature. This form is driven by previous chemical models of CH$_3$OH and CH$_3$CCH, in which prominent abundance enhancement is seen around the two $T_{jump}$ temperatures (see Appendix F). The best-fit parameters $T_{jump}$, $X_{out}$, and $f_{inc}$ for this model (hereafter model A) are listed in Table 7 (Col. A). For all sources we find that with the assumed density profile of $\rho_{bulk}(r)$ the models cannot reproduce the observed high ratios between the higher and lower K components of CH$_3$OH lines, as shown in Fig. 13 (presenting the comparisons between modeled results and observations toward clump G08a and G08b), which points to, as also indicated by the RADEX modeling results, a much higher gas density regime from which these CH$_3$OH higher K components originate. Therefore, we complemented the LIME models with density profiles of the dense gas component $\rho_{dense}(r)$ as in Eq. (4) (the
<figure><img src="image_11.png" /><figcaption>Fig. 9. CH₃OH derived n(H₂) maps from RADEX modeling of all target sources. The beam of the CH₃OH 5-1-4-1 E line is indicated in the bottom left corner. The gray contours indicate the SMA 1.2 mm continuum level from 0.1 to 0.9 × peak flux represented by five steps of uniform interval.</figcaption></figure>

RADEX results of gas density radial profiles):

$$ \rho_{\text{dense}}(r) = \begin{cases} \frac{1}{f_r}\rho_0 & (r < r_0) \\ \frac{1}{f_r}\rho_0(r/r_0)^{q_{\text{radex}}} & (r_0 < r < R_{\text{max}}). \end{cases} \quad (6) $$

Here, r₀ denotes the reference radius of 0.1 pc or 0.05 pc (for G13 and G31) and ρ₀ is the reference density at r₀. These values, together with Rmax, were derived by RADEX modeling. Here fᵣ is a reduction factor applied to ρ₀. This parameter is empirically added to the density profile so that the LIME models more closely match the observed data. In essence, it means that RADEX results of one-component non-LTE modeling tend to overestimate the projected radial density in a 3D structure clump. We manually adjusted fᵣ and Xmol(r) to seek for better fits to the observational data. In what remains, we refer to qradex as qdense as this slope is fitted based on n(H₂) maps of CH₃OH RADEX modeling and is retained as the slope for the dense gas profile in full radiative transfer LIME models. The best-fit model parameters Tjump, Xout, finc, and fᵣ for this model (hereafter model B) are summarized in Table 7 (Col. B). Figure 14 shows a comparison between the CH₃OH line profiles reproduced from model B and the observations toward clump G08a and G08b. The comparisons between model A, model B, and the observations for other target sources are shown in Appendix G.

### 3.7. Comparison of the samples: density and temperature structures

We make a comparison of the fitted and refined T(r) profiles (Eq. (3)) of all sources in Fig. 15 (left panel). We can see (for all the clumps) that at 0.1 pc the resolved gas temperatures are in the range 30–80 K, and at 1 pc around 20–30 K. Moreover, the gas temperature at a certain clump radius is not a monotonic function of the bolometric luminosity of the clump. Particularly, the hot massive core G31 and the source G13 display higher temperatures in the inner regions than their immediate more luminous sources in the sample. We discuss these temperature profiles in more detail in Sect. 4.1.

In Fig. 15 (right panel), we also compare the derived radial gas density profile of the dense gas from CH₃OH modeling with RADEX/LIME. It can be seen that n(H₂) is several ~10⁵ to 10⁷ cm⁻³ at ~0.2–0.3 pc (projected) radii. In the inner ~0.1 pc where the SMA identified continuum cores (Sect. 3.2), n(H₂) ranges from several ~10⁶ to 10⁸ cm⁻³. There are exceptionally high n(H₂) values at the center of G31. Although we have verified the high level of dense gas of this source compared with the rest of the sample by full radiative transfer modeling of CH₃OH lines (Appendix F), we caution that in this density regime the CH₃OH (5-4) lines become heavily optically thick and the critical density for the considered line transitions is reached (e.g.,
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Radial intensity profile comparisons between observations and best-fit RADMC-3D models. Gray horizontal dashed lines indicate the noise level ($3\sigma$). Gray vertical lines indicate the clump radius used in the modeling. The dotted line indicates beam shape in each plot. For source G18, G19, G08a, G13, and G31, the model fit after re-adjusting $T(r)$ is shown. The gray vertical line indicates the clump radius $R_{\text{clump}}$.</figcaption></figure>

$n > n_{\text{crit}}$, Table 3), such that the relative differences between the level populations do not serve as ideal densitometers anymore. Nonetheless, we can safely argue that the hot massive core G31 has much higher gas densities in its inner region than other sources, which is also reflected by the very monolithic nature of its central core from higher angular resolution observations (~2000 au, see Beltrán et al. 2018). For source G10, there are prominently higher gas densities at extended radii of 0.1–0.4 pc than for other sources, which is related to the presence of a large disk-like flattened structure. We discuss further the density profiles among the sample in Sect. 4.2.

We compare the steepness of the radial gas density profiles ($q_{\text{bulk}}$ and $q_{\text{dense}}$, for $\rho_{\text{bulk}}$ and $\rho_{\text{dense}}$) as a function of source

evolutionary stages, which is indicated by the clump bolometric luminosity-to-mass ratio $L/M$ (Fig. 16). Using the Spearman correlation measure, we find that there are positive correlations (correlation coefficient $\rho = -0.95$ and $-0.61$) between the density power-law slopes with $L/M$ for the dense gas component and for the bulk gas structures, although the significance of the correlation of the former is low ($p$-value=0.15). The slopes range from $-0.6$ to $-1.7$ for the bulk gas, and $-0.25$ to $-1.7$ for the dense gas, for $L/M$ spanning from 10 to ~100 ($L_{\odot}/M_{\odot}$) of all sources. The correlation between $L/M$ and the density slope representing the bulk gas distribution is clearly stronger. A similar evolutionary trend was reported by Beuther et al. (2002) based on the analyses of 1.2 mm dust continuum emission of a sample
<figure><img src="image_13.png" /><figcaption>Fig. 11. Comparison of SEDs of the best-fit RADMC-3D models with measured multi-wavelength fluxes (green dots with error bars indicating 0.8 and 1.2 times the flux level) for each source. The black line indicates the SED generated from assumed T(r) and the corresponding best density profile fits. The blue dashed line indicates the SED generated from refined T(r) and the re-iterated best density profile fits. The blue shaded regions indicate a 20% difference around the blue dashed SED profile. The red line shows the SED generated by self-consistently calculating the dust temperature adopting a central heating ZAMS star plus the re-iterated best density profile (for more details, see Appendix E).</figcaption></figure>

of massive clumps, in which the bulk gas density structure was probed. Comparably, other works on the density structures of massive clumps typically derived power-law slopes from -2.25 to -0.75 and peaking at -1.8 to -1.6 (Mueller et al. 2002; Beuther et al. 2002; van der Tak et al. 2000). In the early-stage sources (L/M < 20), the slopes we derived are relatively shallow (> -1.0) for both the bulk gas and dense gas density structures. We note that the slope derived for the dense gas structure of the early-stage source G18 is valid for a confined region of ~0.1 pc (just above the beam size), which merely reflects a pocket of dense gas that is rather compact and remains unresolved. Yet, for other early-stage sources (G28 and G19), the statistics for determining the density slope are adequate, and these two sources exhibit shallow slopes of ~-0.6. We note that different analysis methods of density structure could result in systematic biases in the derived density slopes. In addition, the analyses of dust continuum that were based on the optically thin assumption, instead of relying on full radiative transfer models, suffer from the degeneracy of density and temperature profiles in determining the radial intensity profiles. Moreover, close to the source center, the optically thin assumption for dust emission may also break down. Although a qualitative comparison can be made, a careful gauge between the different analyses conducted are necessary for a stringent comparison between different works. In Sect. 4.5 we discuss the relation between density profiles of massive clumps and statistics on cloud structure in greater detail. We also elaborate on the physical implications by comparing ρ<sub>bulk</sub> and ρ<sub>dense</sub> from the sample.

## 3.8. Molecular linewidths and virial parameter

To understand the dynamic states of the target clumps, we examined how the linewidths and virial parameters vary with clump radii. Part of these analyses were based on the thermometer lines, CH<sub>3</sub>CN, H<sub>2</sub>CS, and CH<sub>3</sub>CCH. They primarily trace the dense

gas close to the centers (0.1–0.4 pc) of the clumps. In addition, we examined the CS and C<sup>3</sup>S (5-4) and H<sup>13</sup>CO<sup>+</sup> (3-2) lines which can trace spatially more extended clump structures due to their lower excitation conditions. We performed single component Gaussian fits to the CS, C<sup>3</sup>S (5-4), and H<sup>13</sup>CO<sup>+</sup> (3-2) line cubes in a pixel-by-pixel manner to obtain the linewidth maps. For the analysis, we trimmed the pixels that have fitting errors of linewidth larger than 2 times the velocity channel widths (ΔV < 2 km s<sup>-1</sup>).

The virial parameter α<sub>vir</sub> characterizes an important aspect of the physical states of the molecular clumps. The ordinary definition of α<sub>vir</sub> (i.e., ignoring magnetic field; see Bertoldi & McKee 1992) is

$$ \alpha_{\text{vir}} = a_1 \frac{2T}{|W|}, \qquad (7) $$

where T = 3/2 M<sub>enc</sub>σ<sub>rms</sub><sup>2</sup> is the kinetic energy, W = -3/5 a<sub>1</sub> GM<sub>enc</sub><sup>2</sup>/R is the gravitational potential energy, M<sub>enc</sub> is the enclosed mass, and a<sub>1</sub> is a geometric factor that accounts for the inhomogeneity of the density distribution (e.g., Bertoldi & McKee 1992; McKee & Holliman 1999). For a spherical clump that has a ∝ r<sup>-q</sup> radial gas density profile, a<sub>1</sub> = (1+q/3)/(1+2q/5). With this definition, a source in energy equipartition (T ~ |W|) has a critical virial parameter of α<sub>cr</sub> = 2a<sub>1</sub>. In a virialized source (2T ~ |W|), it stands that α<sub>vir</sub> = a<sub>1</sub> with α<sub>vir</sub>/α<sub>cr</sub> = 0.5. In the following we refer to the states of α<sub>vir</sub>/α<sub>cr</sub> < 0.5, ~0.5–1, and >1 as sub-virial, virial, and super-virial, respectively.

When deriving α<sub>vir</sub>, it is critical that the tracers observed for the measurement of M<sub>enc</sub> and σ<sub>rms</sub> predominantly emanate from the same gas entity (Traficante et al. 2018). The mass tracer we adopted, which is the dust continuum emission, traces a broad range of gas volume density values distributed over a wide range of radius values. Our selected tracers to indicate linewidths, as the way the temperature profile is measured, show
<table><caption>**Table 5.** Parameters of CH₃OH derived radial density $\rho_{\text{dense}}(r)$ and multi-thermometer derived temperature profiles $T(r)$.</caption><thead><tr><th>Source</th><th>ρ<sub>0.1 pc</sub><br>× 10<sup>6</sup> (cm<sup>-3</sup>)</th><th>Power-law slope q<sub>dense</sub></th><th>R<sub>eff</sub><sup>(a)</sup><br>(pc)</th><th>R<sub>max</sub><sup>(b)</sup><br>(pc)</th><th>T<sub>in</sub><br>(K)</th><th>r<sub>in</sub><br>(pc)</th><th>T<sub>out</sub><br>(K)</th><th>r<sub>out</sub><br>(pc)</th></tr></thead><tbody><tr><td>G18</td><td>2.9(0.2)</td><td>-0.26(0.07)</td><td>0.08</td><td>0.10</td><td>64.7(1.1)</td><td>0.016<sup>(e)</sup></td><td>19.2(0.1)</td><td>0.25(0.05)</td></tr><tr><td>G28</td><td>6.3(0.4)</td><td>-0.83(0.07)</td><td>0.20</td><td>0.30</td><td>96.1(1.2)</td><td>0.02</td><td>21.5(0.1)</td><td>0.80(0.10)</td></tr><tr><td>G19</td><td>5.3(0.3)</td><td>-0.61(0.05)</td><td>0.20</td><td>0.25</td><td>119.5(2.5)</td><td>0.008<sup>(e)</sup></td><td>21.2(0.1)</td><td>0.25(0.03)</td></tr><tr><td>G08a</td><td>10.3(1.2)</td><td>-1.33(0.12)</td><td>0.30</td><td>0.50</td><td>128.6(1.8)</td><td>0.008<sup>(e)</sup></td><td>22.4(0.1)</td><td>1.1(0.14)</td></tr><tr><td>G31<sup>(c)</sup></td><td>8660 (3000)</td><td>-3.22(0.37)</td><td>0.34</td><td>0.42</td><td>400.0<sup>(d)</sup></td><td>0.01<sup>(e)</sup></td><td>18.3(1.4)</td><td>0.27(0.03)</td></tr><tr><td>G08b</td><td>22.9(0.7)</td><td>-1.35(0.02)</td><td>0.34</td><td>0.50</td><td>135.0(1.4)</td><td>0.02</td><td>22.3(0.1)</td><td>0.80(0.08)</td></tr><tr><td>G13</td><td>9.1(1.8)</td><td>-1.67(0.20)</td><td>0.20</td><td>0.40</td><td>229.6(23.1)</td><td>0.004<sup>(e)</sup></td><td>17.0(1.5)</td><td>0.20(0.07)</td></tr><tr><td>G10</td><td>51.6(2.6)</td><td>-1.07(0.03)</td><td>0.28</td><td>0.42</td><td>172.4(4.8)</td><td>0.02</td><td>24.7(0.6)</td><td>0.90(0.31)</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Effective radius $R_{\text{eff}}$ is defined as $\pi R_{\text{eff}}^2 = A$, where the A is the CH₃OH emission area in which $n(H_2)$ can be reliably derived. <sup>(b)</sup>The largest radius (distance to the center) of $n(H_2)$ map derived by CH₃OH; this is due to the irregular shape of the emission area. <sup>(c)</sup>The density profile of source G31 is better described by a two-component power-law form consisting of a shallow slope in the inner region followed by a steep slope of -5.22(0.045) in the outer region. <sup>(d)</sup>Upper limit is set to 400 K in the fit. <sup>(e)</sup>For sources G18, G19, G08a, G31, and G13, $r_{\text{in}}$ is a re-adjusted parameter based on SED calculation and comparison presented in Fig. 11 (temperature profile shown as blue lines in Fig. 20).

<table><caption>**Table 6.** Best-fit parameters from RADMC-3D modeling of the dust continuum in 350 or 450 µm, and 870 µm.</caption><thead><tr><th>Source</th><th>ρ̄<sup>(a)</sup><br>(10<sup>3</sup> cm<sup>-3</sup>)</th><th>q<sup>(b)</sup></th><th>R<sub>clump</sub><br>(pc)</th><th>ρ<sub>0.1 pc</sub><sup>(c)</sup><br>(10<sup>4</sup> cm<sup>-3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>3</sup> M<sub>☉</sub>)</th><th>L<sub>bol</sub><sup>(d)</sup><br>(10<sup>4</sup> L<sub>☉</sub>)</th><th>M<sub>tot</sub> (< 0.5 pc)<sup>(e)</sup><br>(10<sup>3</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>G18</td><td>14.8</td><td>-0.57</td><td>0.8</td><td>3.81</td><td>1.9</td><td>0.8</td><td>0.7</td></tr><tr><td>G28</td><td>12.3</td><td>-0.66</td><td>1.25</td><td>5.1</td><td>6.9</td><td>6.1</td><td>0.8</td></tr><tr><td>G19</td><td>7.8</td><td>-1.37</td><td>1.04</td><td>10.4</td><td>2.5</td><td>3.1</td><td>0.7</td></tr><tr><td>G08a</td><td>13.2</td><td>-1.19</td><td>1.00</td><td>12.3</td><td>3.8</td><td>3.5</td><td>1.1</td></tr><tr><td>G31</td><td>11.7</td><td>-1.22</td><td>2.18</td><td>29.9</td><td>35.3</td><td>57.0</td><td>2.6</td></tr><tr><td>G08b</td><td>3.6</td><td>-1.66</td><td>2.00</td><td>22.9</td><td>8.2</td><td>24.3</td><td>1.3</td></tr><tr><td>G13</td><td>9.5</td><td>-1.41</td><td>0.88</td><td>10.8</td><td>1.9</td><td>3.9</td><td>0.8</td></tr><tr><td>G10</td><td>7.8</td><td>-1.52</td><td>1.58</td><td>23.4</td><td>8.9</td><td>65.3</td><td>1.6</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Average density within clump radius of $R_{\text{clump}}$. <sup>(b)</sup>Density power-law slope. <sup>(c)</sup>Density at 0.1 pc. <sup>(d)</sup>Bolometric luminosity calculated by using the SED profile shown in Fig. 11 (blue lines for source G18, G19, G08a, G31, and G13). <sup>(e)</sup>Total mass within 0.5 pc from clump center.

emission of progressively larger radii, which are complemented with the aforementioned three tracers of more extended emission. We can now examine the spatial variation of linewidths and $\alpha_{\text{vir}}$ based on multiple tracers that cover distinct critical densities (Table 3), and hence different spatial scales. We evaluated how $\alpha_{\text{vir}}/\alpha_{\text{cr}}$ varies with radius using the best-fit density models from the RADMC-3D continuum modeling (Appendix E) to obtain $M_{\text{enc}}$, $a_1$ and the linewidth maps from the tracers. Compared with the SMA observations, we note that the RADMC-3D models constrained by the coarser resolution single-dish continuum data systematically underpredict the 1.2 mm fluxes in the inner radii. To avoid this bias, we adopted the $M_{\text{enc}}$ as $M_{\text{Abel core}}$ (Sect. 3.2) for the inner regions. We recall that $M_{\text{Abel core}}$ is calculated by applying the derived $T(r)$ to SMA 1.2 mm continuum. We discuss the obtained radial profiles of linewidth and virial parameter in Sect. 4.3.

### 3.9. Molecular abundance and abundance ratios

To facilitate the analysis on clump evolutionary stages, we derived the LOS integrated abundance maps ($N_{\text{mol}}/N(H_2)$) for some relevant molecular species for all sources. The bulk gas density profiles ($\rho_{\text{bulk}}(r)$, Sect. 3.6, and Appendix E) were adopted and smoothed to the angular resolution of the specific line transition when deriving $N_{\text{mol}}/N(H_2)$. The calculation

of $N_{\text{mol}}$ for CH₃CCH, CH₃CN, H₂CS, and CH₃OH is introduced in Sects. 3.4 and 3.5. The calculations of the $N_{\text{mol}}$ maps of the CS/C³⁴S, SO, SO₂, and CCH lines were based on LTE assumption and are detailed in Appendix D. We then derived the projected radial averaged abundance profiles for each molecule. Naturally, the projected radial averaging suppresses the contrast in the spatial variations of the abundance for the molecules that are enriched in the clump center or other localized regions (which reduces $N(H_2)$ to localized values rather than integration along the LOS extension). Nevertheless, this does not qualitatively change the overall radial trends as long as the abundance is increasing or decreasing with radius monotonically, with a steeper profile than that of the column density, while the latter is rather shallow, following $\Sigma \propto \rho r \propto \rho^{1+q} < \rho^{-0.7}$. In Figs. 17 and 18, we show the relations between the clump bolometric luminosity and the radial abundance variations of the carbon-chain and sulfur-bearing molecules in consideration, respectively.

In Fig. 17, we note that the abundances of CH₃OH, CH₃CCH, and CH₃CN show similar behaviors. They appear largest in the hot massive cores G31; in the rest of the sources, the abundances of CH₃OH and CH₃CCH are in the range of several $10^{-9}-10^{-8}$, while the abundances of CH₃CN are in the range of $10^{-10}-10^{-9}$. These abundances are slightly positively correlated with the clump luminosity. The abundances of CH₃CN and
<table><caption>**Table 7.** Best-fit CH<sub>3</sub>OH and CH<sub>3</sub>CCH abundance results of LIME modeling based on density model from A: RADMC continuum modeling as listed in Table 6; B: manually adjusted RADEX radial density profile as listed in Table 5.</caption><thead><tr><th rowspan="2">Source</th><th colspan="3">A</th><th colspan="3">CH<sub>3</sub>CCH</th><th colspan="3">B</th><th colspan="3">CH<sub>3</sub>CCH</th></tr><tr><th>X<sub>out</sub><br>(×10<sup>-10</sup>)</th><th>T<sub>jump</sub><br>(K)</th><th>f<sub>inc</sub><br>(<sup>a</sup>)</th><th>X<sub>out</sub><br>(×10<sup>-10</sup>)</th><th>T<sub>jump</sub><br>(K)</th><th>f<sub>inc</sub></th><th>f<sub>r</sub><br>(<sup>b</sup>)</th><th>X<sub>out</sub><br>(×10<sup>-10</sup>)</th><th>T<sub>jump</sub><br>(K)</th><th>f<sub>inc</sub></th><th>X<sub>out</sub><br>(×10<sup>-10</sup>)</th><th>T<sub>jump</sub><br>(K)</th><th>f<sub>inc</sub></th><th>f<sub>f</sub>dens<br>(<sup>c</sup>)</th></tr></thead><tbody><tr><td>G28</td><td>10</td><td>30</td><td>10</td><td>10</td><td>30</td><td>10</td><td>3</td><td>1.5</td><td>80</td><td>150</td><td>0.018</td><td>30</td><td>200</td><td>0.13</td></tr><tr><td>G19</td><td>15</td><td>80</td><td>10</td><td>75</td><td>—</td><td>1</td><td>3</td><td>1.5</td><td>80</td><td>100</td><td>5</td><td>80</td><td>30</td><td>0.19</td></tr><tr><td>G08a</td><td>10</td><td>80</td><td>20</td><td>100</td><td>80</td><td>20</td><td>3</td><td>0.84</td><td>80</td><td>35</td><td>6</td><td>—</td><td>1</td><td>0.06</td></tr><tr><td>G31</td><td>2.5</td><td>80</td><td>10</td><td>20</td><td>80</td><td>10</td><td>50</td><td>4.0</td><td>120</td><td>8</td><td>12</td><td>120</td><td>2</td><td>0.09</td></tr><tr><td>G08b</td><td>10</td><td>80</td><td>5</td><td>100</td><td>—</td><td>1</td><td>5</td><td>0.6</td><td>80</td><td>40</td><td>5</td><td>80</td><td>70</td><td>0.05</td></tr><tr><td>G13</td><td>10</td><td>80</td><td>5</td><td>20</td><td>80</td><td>10</td><td>3</td><td>0.67</td><td>80</td><td>4</td><td>1.9</td><td>80</td><td>3</td><td>0.28</td></tr><tr><td>G10</td><td>10</td><td>80</td><td>5</td><td>—</td><td>—</td><td>—</td><td>7.5</td><td>0.54</td><td>40, 100</td><td>4, 75</td><td>—</td><td>—</td><td>—</td><td>0.05</td></tr></tbody></table>

Notes. “—” denotes an invalid or unavailable parameter. (<sup>a</sup>) *f*<sub>inc</sub> represents increase factor of the abundance jump model (Eq. (5)). (<sup>b</sup>) Reduction factor applied to RADEX density results in the modeling (Eq. (6)). (<sup>c</sup>) Dense gas volume filling factor defined in Appendix F.

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Projected radial averaged *n*(*H*<sub>2</sub>) radial profiles derived from *n*(*H*<sub>2</sub>) maps shown in Fig. 9. The thick gray line indicates the best-fit single power-law model (beam convolution considered). The gray shadowed band indicates the 3σ confidence band of the best-fit model. The model parameters and 1σ errors are listed in Table 5.</figcaption></figure>

CH<sub>3</sub>OH appear more tightly correlated with the source temperatures in the inner ~0.1 pc regions (see the insets in Fig. 17). These trends are consistent with theoretical predictions that the de-sorption of these molecules from the grain surface becomes more efficient with increasing temperature. We note that these results cannot be obtained if the gas temperature distributions are derived based merely on the assumptions of bolometric luminosity scaling instead of being derived based on multi-transition rotational temperature maps, since we have previously seen that the gas temperatures at certain radii do not necessarily increase monotonically with clump luminosity (Fig. 15, left panel). This result demonstrates the importance of measuring detailed temperature profiles when studying chemical evolution.

We observed a weak correlation of abundance of CH<sub>3</sub>CCH and bulk gas temperature, which is reminiscent of the small variation of the CH<sub>3</sub>CCH abundance toward massive clumps of various evolutionary stages reported by Giannetti et al. (2017) (see also Öberg et al. 2014). Other higher angular resolution observations of this species indicated a mixed behavior of its spatial distribution, depending on whether the emission coincides with the localized hot cores or appears offset and/or shows more extended structures (Bøgelund et al. 2019, Öberg et al. 2014, Fayolle et al. 2015). Comparing rotational temperature maps of CH<sub>3</sub>CCH (12-11), the temperatures of G13 and G31 in the core region is 20–30 K higher than those of G19, G08a, and G08b. In addition, the higher J transitions CH<sub>3</sub>CCH (14-13), (15-14) trace systematically higher temperatures in the core region toward all sources. This evidence indicates that the emission of CH<sub>3</sub>CCH has a contribution from gas components residing inside or in the vicinity of hot cores. As for the result of CCH, except for the earliest-stage source G18, its abundance appears enhanced at outer radii. We can also see that the abundance of CCH measured in the inner ~0.1–0.2 pc region is anti-correlated with the clump luminosity. These results are consistent with previous observations that show shell-like CCH emission toward late-stage massive star-forming regions (Beuther et al. 2008; Jiang et al. 2015).

The abundance ratios of CCH, CH<sub>3</sub>CCH, and CH<sub>3</sub>CN with CH<sub>3</sub>OH in the clump center (0.1–0.15 pc) are shown in Fig. 19. Here, normalization with CH<sub>3</sub>OH allows exploring the chemical evolution or initial condition by eliminating the effect of different desorption levels among the sample. There are substantial variations of [CCH]/[CH<sub>3</sub>OH] and [CH<sub>3</sub>CN/CH<sub>3</sub>OH] across *L/M*, while variations of [CH<sub>3</sub>CCH/CH<sub>3</sub>OH] appear moderate.
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** LIME modeling result (best-fit parameters listed in Table 7, Col. A) based on best-fit density model from RADMC-3D continuum modeling. From left to right: Annular beam-averaged spectra from the continuum center to the outer envelope. Considering the typical beam FWHM of our observations: the distance from the center of each annular region to the center of the source is indicated at the top of each spectra. The line components of A- and E-type CH₃OH are indicated with short-dashed vertical lines in green and gray, respectively.</figcaption></figure>

The abundances of sulfur-bearing species in the central regions of the clumps also show correlation with temperatures (Fig. 18). Comparing the relative abundance of C³⁴S and H₂CS with respect to SO as a function of *L/M* (Fig. 19, bottom panel), it seems there is a mixed behavior; except for source G13, the other sources show a weak increasing trend above *L/M* ≥ 10. For *X(SO₂)/X(SO)*, there is a moderate increasing trend with source *L/M*. We discuss the abundance variations further in Sect. 4.4 in a broader context, with comparisons with published results from chemical modeling.

# 4. Discussion

## 4.1. Temperature structure and heating mechanisms of massive star-forming clumps

The temperature measurement from multiple $T_{\text{rot}}$ maps and the fitted and refined temperature profile $T(r)$ (Eq. (3)) are shown

in Fig. 20. The dlog*T*/dlog*R* profiles are also summarized in the bottom right panel. Except for sources G13 and G31, the dlog*T*/dlog*R* profiles asymptotes from −0.5 to zero from inner to outer radii. The levelling-off of the temperature at the outer radius of the clump is expected since all these massive star-forming clumps are immersed in intense interstellar radiation fields. At gas densities >10⁴.⁵–10⁵ cm⁻³, thermal coupling between gas and dust can be quickly achieved due to collisions (Goldsmith 2001; Glover & Clark 2012). However, the 10⁴.⁵ cm⁻³ density threshold is not met in the outer layers of the sources that have lower masses. As we explain below, the thermal decoupling between gas and dust is seen in some of these sources. The outer envelope dust temperature for all these sources flattens around 18–25 K. These values are consistent with an elevated infrared radiation field associated with these regions. Based on these temperatures, the scaling factors of the local ISRF (Mathis et al. 1983) to characterize the radiation field surrounding these
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** LIME modeling result (best-fit parameters listed in Table 7, Col. B) after manually adjusting the density profile obtained from RADEX modeling (Eq. (4)), which is prescribed as a piecewise power law (Eq. (6)). From left to right: annular beam-averaged spectra from the continuum center to the outer envelope. Considering the typical beam FWHM of our observations: the distance from the center of each annular region to the center of the source is indicated at the top of each spectra.</figcaption></figure>

clumps are roughly $\ge 10^{2.5}-10^4$ (estimated in the optically thick limit, Krumholz 2014), which are typical values found in the vicinity of Galactic massive star-forming complexes (Binder & Povich 2018).

The $T(r)$ and $d\log T/d\log R$ in the inner regions can be approximated by the analytic temperature profile of thermally balanced dust grains distributed around a central heating point source (e.g., Adams 1991). With optically thin condition and a submillimeter dust opacity spectral index ($\beta$), the model of Adams (1991) predicted that the radial temperature profile for dust grains in thermal balance around a central heating point source is proportional to $r^{-2/(4+\beta)}$. If $\beta$ has no spatial variations, then $d\log T/d\log R$ should be a constant of radius. In the diffuse interstellar medium $\beta$ is around 1.8 (for a review, see Hildebrand 1983), which yields a temperature slope of $\sim -0.35$. In high-density regions, $\beta$ may become lower due to dust growth (e.g.,

Ossenkopf & Henning 1994), resulting in a steeper temperature profile. Values of $\beta$ lower than 1.8 have been reported by some previous observations toward Class 0-II young stellar objects (YSOs), and toward protostellar and prestellar cores (e.g., Beckwith & Sargent 1991; Jørgensen et al. 2007; Bracco et al. 2017, and references therein). This comes with the caveat that the previous (sub)millimeter observations of dust growth might have systematically underestimated $\beta$ values owing to underestimating optical depths (e.g., Li et al. 2017), neglecting the effect of dust scattering opacity (e.g., Liu 2019), and the effects of temperature mixing when performing SED fittings (e.g., Juvela et al. 2018). From a modern point of view, on the spatial scales of molecular clumps, there might not yet be solid evidence for the presence of $\beta < 1$ (i.e., a flattened SED curve at longer wavelengths). We note that observations revealing a prevalent excess of 3 mm emission compared to the generic modified blackbody model have been
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Radial temperature and density profiles of target clumps. **Left panel:** radial temperature profiles $T(r)$ for all the target sources (Eq. (3)), refined $T(r)$ is used for the relevant source). The thickness of the lines increases with increasing luminosity. **Right panel:** radial (projected radial averaged) density profile of the dense gas ($\rho_{\text{dense}}$) of all the target sources, from $n(\text{H}_2)$ maps derived by RADEX modeling of CH$_3$OH lines. The range of the y-axis is trimmed to increase contrast. For both panels the luminosity and clump mass are calculated from RADMC-3D best-fit model (Table 6), which takes into account all the gas components present in the clumps.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Radial gas density slopes of all target clumps. **Left panel:** density power-law slope derived from continuum ($q$) based on RADMC-3D modeling detailed in Appendix E. The luminosity, clump mass, and enclosed mass within 0.5 pc are calculated from the RADMC-3D best-fit model, as listed in Table 6. **Right panel:** density power-law slope derived from $n(\text{H}_2)$ maps ($q_{\text{radex}}$) from the CH$_3$OH RADEX modeling detailed in Sect. 3.5. The yellow horizontal line in both plots shows a slope of -1.5, indicating the free-falling density profile of a singular isothermal sphere (with an initial density slope of -2) as in Shu 1977, and the attractor solution of the gravoturbulent collapsing in Murray et al. (2017).</figcaption></figure>

reported (e.g., Lowe et al. 2021); however, the origin of such a feature remains uncertain.

In Fig. 20, we also compare the observed temperature profiles with the centrally heated models evaluated for $\beta = 1$ (i.e., $T = 2.70(L/L_\odot)^{0.2} (R/1\text{pc})^{-0.4}$), based on the assumption of a dust sublimation temperature of $1.1 \times 10^3$ K in the derivation (Adams 1991, for more calculation details see the Appendix therein), and with the simple expectation based on the Stefan-Boltzmann law (i.e., $T = 0.86(L/L_\odot)^{0.25} (R/1\text{pc})^{-0.5}$). We found that these profiles qualitatively agree with the measurements of $T(r)$ except for G13 and G31.

In Fig. 20, it can be seen that the temperature profiles of G13 and G31 deviate from the form of a single power law. Specifically, they show abrupt changes as well as elevated temperatures

in the ranges of 0.1–0.3 pc and 0.1–0.5 pc radii, respectively. After we adjusted $T(r)$ (Appendix E) according to dust SED profiles, G13 shows a less prominent temperature enhancement in the center, while G31 still stands out. Beltrán et al. (2018) also noted the steep temperature profile of G31 within the central 0.1 pc. Their measurement consistently falls onto the functional form we fitted (i.e., the thick purple line in Fig. 20) which was based on our independent measurements of $T_{\text{rot}}$ on the larger spatial scales.

The rapid decrease in radial temperatures observed in G31 and G13 may be explained by their density profiles of the embedded dense gas structures (i.e., $\rho_{\text{dense}}(r)$). In Fig. 12, it can be seen that these two sources have the most steeply decreasing $\rho_{\text{dense}}(r)$. In addition, their central ~0.1 pc regions show prominent
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Abundance profiles of molecules CH<sub>3</sub>OH, CH<sub>3</sub>CN, CCH, and CH<sub>3</sub>CCH toward target sources. The symbols are color-coded based on relative distance to the clump center (the 1.2 mm continuum peak). Stars represent abundances at the clump center, and hollow circles or triangles indicate abundances obtained from different radii: the larger and darker these symbols, the closer the distance to the continuum peak. For clump G28 the data points at outer radii are shown as triangles, whereas those for other clumps are shown as circles, to further avoid confusion. X(CCH) for source G10 is taken from Jiang et al. (2015). The inset plot for CH<sub>3</sub>CN shows the central abundance vs. gas temperature at 0.1 pc; the plot for CH<sub>3</sub>CCH shows the envelope abundance vs. gas temperature at 0.3 pc.</figcaption></figure>

high-density plateaus. The optically thin assumption of Adams (1991) may break down in the central region of the G31 and G13. This higher concentration of dense gas may steepen the temperature distribution in the inner envelope according to $\propto r^{-1-q/4-\beta}$ (e.g., Adams & Shu 1985; Rolffs et al. 2011), where $q$ is the power-law index of radial density profile. This is close to the diffusion approximation, which effectively means that a higher optical depth gas would mimic the lower value of $\beta$ in determining the temperature structure. The presence of flattened (protostellar) disks could also induce a steeper gradient (~-0.75) of temperature variations due to the gas heating by infall and accretion shocks (e.g., Lynden-Bell & Pringle 1974; Walch et al. 2009). Interestingly, hydrodynamic calculations of protostellar collapse demonstrate that a centrally flattened density profile results in a transitory phase of energetic accretion (Foster & Chevalier 1993, see also Henriksen et al. 1997), which may also be related to the elevated temperatures of G13 and G31.

Observations toward dense massive cores and envelopes of YSOs and disks generally find temperature slopes in the range of ~[-0.35, -0.7] (e.g., Palau et al. 2014; Beuther et al. 2007; Persson et al. 2016; Jacobsen et al. 2018; van't Hoff et al. 2020; Gieser et al. 2021). However, the exceptionally steep temperature profile of slope steeper than -0.9 is seen at 1000-2000 au around the massive YSO W3IRS4 (Mottram et al. 2020), reminiscent of the result of G31 based on observations of similar spatial scales (Beltrán et al. 2018).

Strong radiative feedback has been invoked as a possible mechanism to prevent overfragmentation, which favors the formation of massive stars and such effect can be influential in ultra-dense environment (e.g., Krumholz et al. 2012). When observed with ~2000 au resolution, source G31 consists of two cores with one major core dominating the emission (~60 times flux difference at 1 mm continuum, Beltrán et al. 2018). The elevated temperature in the inner 0.1 pc might have suppressed the fragmentation of the envelope structure of the main core, leaving only one satellite core surrounding it. The highly concentrated dense gas structure of G13 traced by methanol lines (Fig. 9) that does not extend further beyond its continuum emission (except to the west) may also indicate a featureless fragmentation; high-resolution (1") mid-IR imaging by Varricatt et al. (2018) reveals a binary system, although the mass contrast between the two YSOs embedded is not as drastic as that of G31.

## 4.2. Density structure evolution: comparison with theoretical predictions

Figure 9 shows the $n(H_2)$ maps that were derived based on the RADEX modeling for the CH<sub>3</sub>OH lines (Sect. 3.5). In general, $n(H_2)$ decreases radially with respect to the continuum center, although some sources appear notably asymmetric and present localized overdensities at large radii. We use single power-law forms (Eq. (4)) with the parameters listed in Table 5 to describe
<figure><img src="image_20.png" /><figcaption>Fig. 18. Same as Fig. 17, but for H₂CS, C³⁴S, SO, and SO₂.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19. Comparisons of abundance ratios for relevant molecular species of target clumps. Top: abundance ratios of carbon-chain molecules; Bottom: abundance ratios of sulfur-bearing molecules at clump center (0.1–0.15 pc, beam averaged), shown against source L/M ratios.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Derived radial averaged temperature profiles of the target sources from multiple thermometers. The error bars show the standard deviations for each annular average. The dashed and dotted lines show a radial temperature profile that follows $\propto L^{0.25}r^{-0.5}$ ($\beta = 0$) and $L^{0.2}r^{-0.4}$ ($\beta = 1$), respectively (for details, see Sect. 3.6). For G13 and G31 an additional radial temperature profile of $\propto L^{0.17}r^{-0.34}$ ($\beta = 1.8$, Adams 1991) is shown (dash-dotted line). Whenever available, temperature measurements from higher angular resolution observations from previous works are included in the plots as gray crosses. The thick purple line indicates the fitted temperature profile $T(r)$ described in Sect. 3.6. The blue thin line indicates the refined temperature profile by varying $r_{\text{in}}$ in $T(r)$ (Eq. (3)) to fit with dust SED (refined $T(r)$, Appendix E). The plots in the bottom right panel show the first derivative $d \log T/d \log R$ (for all sources, top) calculated from the fitted profile $T(r)$, and the refined $T(r)$ (for five sources, bottom).</figcaption></figure>

the radial change of $n(H_2)$ maps. Comparing the observed profile to the fitted single power-law form (Fig. 12), it can be seen that there are higher density plateaus at the centers of G13 and G31, such that these two sources are better described by piecewise power law including a central density profile of slope ~0. Source G10 harbors even more dense gas at >0.1 pc radii compared to its best-fit single power-law model, which is due to the highly flattened dense gas geometry of an edge-on rotational disk, as revealed by gas kinematics (Liu 2017). For G10 the slope of the dense gas profile $q_{\text{dense}}$ is an average of two distributions: a geometrically flattened high-density pseudo-disk with a slope shallower than -0.5, and an envelope whose gas density decreases sharply (Fig. 9). Clumps G28, G19, and G08a may also harbor high-density plateaus at the centers although they are not as significantly resolved as those in G10, G13, and G31.

We note that the two sets of power-law density slopes derived reflect clump gas in different density regimes: from the single-dish continuum the slope represents density distribution of averaged (mass-weighted) bulk gas ($\bar{\rho} \sim 10^{4.5} \text{ cm}^{-3}$) where the spherical symmetry might be a good approximation of source geometry, while the slope deduced from the methanol emission

reflects the truncated central dense gas portion ($\rho \gtrsim 10^{6.5} \text{ cm}^{-3}$) with an extension of 0.1–0.5 pc. This dense gas structure is likely highly fragmented and of reduced dimension due to the dominant role of self-gravity at increasingly smaller scales and denser regimes.

Figure 16 reveals the evolutionary trend of the radial density profiles of bulk gas and dense gas structures of the target massive clumps. Our results indicate that in the initial stage of massive clump evolution, gas may be less concentrated than the singular isothermal sphere (SIS, Shu 1977), and even the logatropic model (McLaughlin & Pudritz 1997). Observations toward low-mass cores at early stages have generally found density profiles with a central plateau (e.g., Ward-Thompson et al. 1999; Bacmann et al. 2002). This density structure implies low-pressure gradients at small radii, and hence in quasi-static models additional support of magnetic fields is required. The magnetic field and subsequent ambipolar diffusion may explain the existence and evolution of the shallower radial density profile of the gas, which is at an initial stage of equilibrium (e.g., Mouschovias & Morton 1991). Alternatively, in an isothermal state, the incoherently converging compression wave of shocks
<figure><img src="image_23.png" /><figcaption>**Fig. 21.** Derived averaged radial linewidth and virial parameter profiles. The dash-dotted line in each upper panel indicates the power-law fit to the linewidths from the thermometer lines. In each upper panel the gray dotted line indicates the relation found by Caselli & Myers (1995) of Orion low-mass cores scaled up by a factor of 5, which roughly matches the observed linewidth in our case, of $\Delta v/1 \text{ km s}^{-1} = 7.5(r/1 \text{ pc})^{0.53}$. The gray dashed line shows a relation of $\Delta v/1 \text{ km s}^{-1} = 6.0 (r/1 \text{ pc})^{0.2}$. These two reference lines are identical in all plots. In each lower panel, the ratio of $\alpha_{\text{vir}}$ to the critical virial parameter $\alpha_{\text{cr}}$ is shown; the horizontal solid and dotted line indicate a ratio of 1 (equipartition) and 0.5 (virial equilibrium), respectively. All plots except for G10 and G18 share the same legend shown in the middle panel of the last row, where the color-coding for the thermometer lines is the same as in Fig. 20.</figcaption></figure>

propagating outside-in can also disrupt the centrally peaked density profile (Whitworth et al. 1996). Finally, hydrostatic models of a self-gravitating gas externally heated (e.g., Falgarone & Puget 1985) can also result in a similar density configuration. On the other hand, we cannot rule out the possibility that the shallower density profiles reflect the underlying fragmentation that is distributed widely within the clump, and that the fragments do not show significant mass segregation (Sanhueza et al. 2019). Their imprints would result in a uniform density profile for the clumps.

Following the definitions in Sect. 3.6, we denote the power-law slopes for the radial gas density of bulk gas and dense gas as $q_{\text{bulk}}$ and $q_{\text{dense}}$. In G18 and G28, both $q_{\text{bulk}}$ and $q_{\text{dense}}$ are shallower than $-1$ (Fig. 15). This may be due to a combination of the heated gas profiles and a higher level of turbulence,

which effectively slow down gravitational collapse at early stages (Fig. 21; more in Sect. 4.3). Although these two sources have $L/M \lesssim 10$, and are classified to be younger than or just reaching the zero age main sequence (ZAMS) phase (Molinari et al. 2008; Giannetti et al. 2017), it is possible that the enhanced energy release of accretion associated with the massive star formation already introduces an observable temperature gradient (Fig. 20).

The more evolved source G08a has a slightly shallower overall distribution of gas ($q_{\text{bulk}} \sim -1.2$) than its dense gas component ($q_{\text{dense}} \sim -1.3$), although the difference is within the errors. Interestingly, source G19 displays an opposite relation between the two slopes ($q_{\text{bulk}} \sim -1.4$, $q_{\text{dense}} \sim -0.6$; Fig. 16). Comparing the virial state of these two sources (Fig. 21), it seems G19 is close to a global state of energy balance, while G08a has excessive kinetic energy in the central 0.2 pc, possibly caused by energy
transfer and induced motions from gravitational collapse (see Sect. 4.3 for further details). If a steeper density profile corresponds to a more dominant role of gravitational collapse, then the kinematic differences inferred from linewidths are compatible with the relation between density slopes on local (dense gas regime) and global (bulk gas regime) scales for these two sources. We discuss this in more detail in Sect. 4.5.

The overall trend seen in the right panel of Fig. 15 can be qualitatively explained by the growing role of self-gravity in the evolution of gas dynamics in massive clumps (Lee et al. 2015; Gómez et al. 2021). This is illustrated in hydrodynamic simulations of star-forming clouds with continuously driven turbulence (Lee et al. 2015); the simulations show that self-gravity plays the dominant role for all changes that happen with density structures on scales up to >0.1 pc. Their results reveal that the power-law slopes around density peaks change from -1.3 to -1.55 before and after the onset of star formation, of 0.75 and 1.25 times global free-fall timescale ($t_{ff}$). By separating the scales based on the dominant mass contributor of gas or forming stars, Murray et al. (2017) demonstrate that close to the (proto)stars the gas density profile has a slope of -1.5 (an attractor solution, see also Coughlin 2017) and in the outer envelope the slope ranges from -1.6 to -1.8. These values are in quantitative agreement with the evolved sources in our sample. The three sources having bulk gas density slopes shallower than -1 (left panel of Fig. 16) may reflect an initial condition close to a constant-density core when gravitational collapse has not significantly altered the gas density profile, for both the bulk gas and dense gas structures (G18 and G28) or only the dense gas component (G19). The tendency of an initial flat inner density profile remaining flat over one $t_{ff}$ determined by the gas central density is also discussed in Henriksen et al. (1997), which is attributed to the fact that the cloud center has the fastest growing velocity mode. The timescale is fleetingly short without further support of turbulence.

### 4.3. Kinematic state of clumps: radial profiles of molecular linewidth and virial parameter

Based on the calculations in Sect. 3.8, Fig. 21 shows the radial profiles of the molecular linewidths ($\Delta V = 2.355\sigma_v$) for individual species. We fitted a power law to the radial linewidth profiles traced by thermometer lines (i.e., excluding the data points from C$^{34}$S and H$^{13}$CO$^+$). In G08a, G31, and G13, the linewidths traced by these lines decrease radially with power-law indices of -0.40, -0.31, -0.56, respectively. Except for the sources G18 and G19, the linewidths on the extended regions traced by C$^{34}$S and H$^{13}$CO$^+$ are systematically larger than those traced by the thermometer lines. Again, for clump G18 we had to rely on pointed observations from the IRAM 30m telescope to tentatively determine the radii of the linewidths obtained from thermometer lines, based on a fixed temperature profile (Sect. 3.4). We therefore note that the analysis for this clump is not based on the same spatially resolved measurement as the other clumps.

In general, the variation in linewidths does not seem to follow a monotonic radial change, as opposed to the $\sigma_v \propto r^{1/2}$ universal law observed in giant molecular clouds and low-mass cores (i.e., Larson's first law; Larson 1981). This is in agreement with the results presented in Izquierdo et al. (2021, Fig. 7), where they show that the linewidth–size relationship is not uniform, but instead depends on the analyzed spatial scale and the physical processes taking place there (see also Hacar et al. 2016).

In high-mass star-forming regions, the observed linewidths are systematically larger, and $\sigma_v$ may have a shallower

power-law relation with scales (e.g., Plume et al. 1997; Caselli & Myers 1995). In addition, the virial equilibrium naturally infers that $\sigma_v$ has a dependence on gas surface density, such that $\sigma_v \propto (\Sigma R)^{0.5}$ (Heyer et al. 2009).

Moreover, recent 1D simulations incorporating gravitation-driven turbulence (i.e., adiabatic heating; see Robertson & Goldreich 2012) found distinct relations between $\sigma_v$ and spatial scales for star-forming clumps, due to the change of the dominant turbulent driving mechanism from the inner to the outer regions (Murray & Chang 2015). It has been found that $v_T \propto r^{-0.5}$ holds within the sphere enclosed by the stellar influence radius, which is defined as the radius where the enclosed gas mass is comparable (1–3 times) to the stellar mass (Murray & Chang 2015; Murray et al. 2017; see also Coughlin 2017; Xu & Lazarian 2020 for the analytic derivations).

As a rough estimate of the stellar influence radius, assuming the $L_{bol}$ of each source is contributed solely by the luminosity of a single ZAMS star, adopting the stellar evolution model of Choi et al. (2016), the stellar mass $M_*$ is estimated to be 12 $M_\odot$ for G19, 15–16 $M_\odot$ for G08a, G13 and G28, and 30 $M_\odot$ for G08b. Grave & Kumar (2009) fit radiative transfer accretion models of YSOs to source G19 and G13 by building SEDs from the near-infrared to the submm, and their derived stellar masses are consistent with our rough estimates. For G10, previous observations suggest a stellar mass of ~200$M_\odot$ (Liu et al. 2010). We note that inferring stellar mass from luminosity assuming the sole contribution from a single star leads to a lower limit, as these clumps are forming a cluster of stars following the IMF. With the same total bolometric luminosity, this corresponds to a higher total stellar mass. In comparison with the central core masses we derived from 1.2mm continuum (Table 4), it seems we marginally resolved the influence radius (about half beam FWHM) for sources G13, G08a, and G10, but a bit further in the case of G08b. This may partially explain why in the inner region of G08b the observed linewidths do not decrease rapidly with radius. In G10 the observed linewidths have shallower variations with radius. This may be due to the presence of the ~0.2 pc scale edge-on rotational disk (Liu 2017), toward which gas settles into coherent rather than random motions.

For the early-stage clumps G18, G19, and G28, the bolometric luminosity is likely dominated by accretion energy. This leads to an overestimation of the stellar mass, such that in reality there appears to be a more drastic difference between the embedded stellar mass and the gas mass at the scale probed, with respect to our rough estimates. Therefore, it is likely that we do not see the decreasing trend of linewidths with radius in their inner regions due to the too coarse resolution to resolve the stellar influence radius. In addition, assuming turbulent velocity scales as the local infall velocity, Coughlin (2017) demonstrates that there is temporal steepening of the radial profile of the rms velocities, due to gravitational field starting to dominate the dynamics of inflowing gas, which changes from $v_T \propto r^{-0.2}$ to $v_T \propto r^{-0.5}$. Hence, it may also be that early-stage clumps have a shallower inner slope that tend to flatten out, further smeared in coarser resolution measurements.

Exterior to the stellar influence radius, simulations found that $v_T \propto r^{0.2}$ (see also Lee et al. 2015). Compared to classical subsonic turbulence following the Kolmogorov law ($v_T \propto r^{1/3}$) or supersonic turbulence ($v_T \propto r^{1/2}$), the shallower scaling relation between linewidth and spatial scale of the region beyond stellar influence radius could be due to additional energy converted partly from gravitational collapse (Ballesteros-Paredes et al. 2011; Xu & Lazarian 2020) and/or kinetic energy of the extended inflow gas transported by radial motions
<table><caption>**Table 8.** Properties related to radial linewidth and virial parameter profiles of the target clumps.</caption><thead><tr><th>Source</th><th>M<sub>★</sub><sup>(a)</sup><br>(M<sub>☉</sub>)</th><th>Slope<sup>(b)</sup></th><th>ΔV<sub>0.1 pc</sub><sup>(b)</sup><br>(km s<sup>-1</sup>)</th><th>Radial variations of α<sub>vir</sub>/α<sub>cr</sub></th></tr></thead><tbody><tr><td>G18<sup>(c)</sup></td><td>11</td><td>—</td><td>—</td><td>>10 in center (&lt;0.1 pc) and ~1–2 at outer envelope (&gt;0.1 pc)</td></tr><tr><td>G28</td><td>23</td><td>-0.01(0.48)</td><td>3.6(0.6)</td><td>fluctuate at ~0.5–5 at all radii</td></tr><tr><td>G19</td><td>18</td><td>-0.35(0.40)</td><td>4.2(0.5)</td><td>fluctuate at ~0.5–5 at all radii</td></tr><tr><td>G08a</td><td>17</td><td>-0.46(0.08)</td><td>4.6(0.3)</td><td>>2 in center (&lt;0.2 pc) and ~0.5–2 at outer envelope (&gt;0.2 pc)</td></tr><tr><td>G31</td><td>64</td><td>-0.30(0.14)</td><td>6.3(0.4)</td><td>~1 in center (&lt;0.1 pc) and ~0.5–1 at outer envelope (&gt;0.1 pc)</td></tr><tr><td>G08b</td><td>26</td><td>-0.09(0.03)</td><td>3.3(0.1)</td><td>≤1 at all radii</td></tr><tr><td>G13</td><td>19</td><td>-0.56(0.13)</td><td>3.7(0.4)</td><td>>2 in center (&lt;0.1 pc) and ~1–2 at outer envelope (&gt;0.2 pc)</td></tr><tr><td>G10</td><td>68</td><td>-0.01(0.13)</td><td>4.7(0.4)</td><td>~0.5 in center (&lt;0.1 pc) and ~1–2 at outer envelope (&gt;0.1 pc)</td></tr></tbody></table>

**Notes.** (a)Estimated (single) stellar mass based on interpolating the ZAMS models of Choi et al. (2016), from the clump $L_{bol}$ listed in Table 6. (b)Fitted slope and linewidth amplitude at 0.1 pc (and errors) of the linewidth radial profile in the inner region. The fitted relation is indicated as the dash-dotted line in the upper panel of Fig. 21. (c)For G18, the parameters cannot be properly constrained in the same manner as the other clumps (see Sect. 4.3) and are omitted.

(Padoan et al. 2020). Klessen & Hennebelle (2010) demonstrate that the conversion efficiency to internal turbulence depends on the density contrast between the accreting entity and the inflowing gas. These may explain why in the outer regions of G08a, G08b, G13, and G10 the observed linewidths increase with radius, having a power-law slope of ≥0.2–0.3.

In light of the radial profiles of $\alpha_{vir}$, it appears that our target clumps can be classified into three types, based on their $\alpha_{vir}/\alpha_{cr}$ ratios: sources G08b, G31, and G10 are in an overall virial to sub-virial state; G18, G08a, and G13 show super-virial states in the center and sub-virial states in the outer envelope; G28 and G19 have uniform $\alpha_{vir}$ fluctuating mostly above $\alpha_{cr}$ (see the detailed description in Table 8). To explain the observed super-virial states in the central regions of some of these clumps, it is required to introduce other line-broadening mechanisms.

The three sources that have large virial parameters in the centers have the lowest central core masses and densities (Table 4). If central cores are treated as a decoupled structure from the clump envelope, the more massive and denser cores being more sub-virial is consistent with previous observations toward massive star-forming regions (e.g., Kauffmann et al. 2013; Liu et al. 2015; Traficante et al. 2018). However, due to the insufficient angular resolution, the estimate of $\alpha_{vir}$ in the clump center may suffer from larger uncertainties, for example there is possibly higher level of clumpiness with complex velocity structures that confuse the linewidth measurement. In addition, the $\alpha_{vir}$ in the centermost is derived from linewidth of CH$_3$CN. CH$_3$CN has been shown to trace the interface of convergent flows in high-mass clumps (Csengeri et al. 2011), which is also enhanced in shocks (Bell et al. 2014) and associated with hot accretion flows (Liu et al. 2015).

Camacho et al. (2016) and Ballesteros-Paredes et al. (2018) have shown that the observed high $\alpha_{vir}$ may not be indicative of pressure confining or gas dispersal. Alternatively, for low column density clumps this picture might be due to highly dynamic externally driven gas accumulation. For high column density clumps it might be a result of neglecting the stellar mass or accreting materials outside the cores in contributing to the gravitational energy. The latter may qualitatively explain the super-virial state in the center of the three sources G18, G08a, and G13, which does not necessarily indicate a halt of collapse. On the other hand, the inclusion of infall velocities to the observed linewidth may also cause the source which is undergoing collapse to appear seemingly unbound (e.g., Smith et al. 2009). Infall motions of dense cores that have large $\alpha_{vir}$

(~5–8) are revealed by Cesaroni et al. (2019), inside a massive cluster-forming clump.

Giannetti et al. (2017) suggest that CH$_3$CCH is a reliable tracer for the kinematics of bulk gas structures in the massive clumps. From the $\alpha_{vir}$ derived based on the CH$_3$CCH lines, it appears that all clumps are in sub-virial or virial status, but at smaller scales in the central region of the clumps the gas can appear super-virial. In addition to the three clumps showing large $\alpha_{vir}$ in the center, other clumps show an evolution of globally (for all radii) decreasing $\alpha_{vir}$ with increasing $L/M$. Ballesteros-Paredes et al. (2018) and Camacho et al. (2020) investigated the evolution of $\alpha_{vir}$ in molecular clouds based on numerical simulations. They showed that the gas structures may be initially assembled due to the large-scale turbulence. Afterward, the assembled gas structures evolve from super-virial to sub-virial status, and finally approach energy equipartition, or further become super-virial after one $t_{ff}$ due to gas expulsion. Our observations do not show such non-monotonic evolution of $\alpha_{vir}$ with $L/M$. This is likely because all target clumps remain in early to intermediate evolutionary stages, which are prior to or in the stage of active star formation when the gas mass dominates the mass budget, preceding the gas dispersal stage.

### 4.4. Comparison with chemical models: carbon-chain molecules and sulfur-bearing species

The formation pathways of molecules and chemical rates are strongly influenced by gas temperature and volume density. In Sect. 3.9, we used the molecular column densities of CH$_3$CCH, H$_2$CS, CH$_3$CN, CH$_3$OH, C$^{34}$S, CCH, SO, and SO$_2$ to derive the molecular abundance distributions and abundance ratios between relevant species. We discuss comparisons between these results with predictions from chemical models in this section.

#### 4.4.1. CH$_3$OH, CH$_3$CCH, CH$_3$CN, and CCH

CH$_3$CCH is designated as a “cold molecule” (Bisschop et al. 2007) since it presents low excitation temperatures compared to other complex organic molecules (COMs) such as CH$_3$CN and CH$_3$OCH$_3$. The formation pathways of CH$_3$CCH include cold gas phase reactions and grain-surface chemistry (Calcutt et al. 2019). The reactions in gas phase involve thermal desorption which requires low temperatures. Calcutt et al. (2019) suggested that CH$_3$CCH is a “gateway” molecule that traces the interface between hot cores-corinos and the colder envelope, and extends further out in the envelope. The inset in Fig. 17 shows measured
CH₃CCH abundances averaged from the outer region (≥0.1 pc) of the clumps with the gas temperature measured at 0.3 pc. These properties represent the bulk luke-warm gas. The small variation of CH₃CCH abundance across our sample indicates that it is not a sensitive molecule to trace the ≥100 K thermal desorption of hot core regions. Its emission instead traces the more extended gas and is a good kinematic tracer for >0.1 pc gas structures within massive clumps (e.g., Giannetti et al. 2017). There is also no drastic radial change of the abundances of the other two COMs CH₃OH and CH₃CN for individual sources (Fig. 17). This is expected since our angular resolution is not sufficiently high to reveal the >100 K region (confined within 0.1 pc) where prominent abundance enhancement of COMs is expected (Garrod et al. 2017).

The abundance profiles of CCH exhibit an increase toward the outer clump envelope for most the sources in the sample (Fig. 17) and an overall decrease in abundance with clump luminosity. These results reflect that CCH is abundant at an early stage and then transformed to other molecules in the clump center (Beuther et al. 2002). The abundance of CCH in outer regions is maintained by reproduction of CO dissociated by an external UV field. Similarly, in early-stage low-mass cores, CCH is also tracing more quiescent gas, e.g., the outer part of a circumbinary envelope (van Dishoeck et al. 1995).

We compare the abundance ratios between CCH, CH₃CN, and CH₃CCH with CH₃OH in Fig. 19. In low-mass star-forming cores, the ratio [CCH]/[CH₃OH] is usually used to distinguish less evolved warm carbon-chain cores with hot corinos (e.g., Graninger et al. 2016). The strong anti-correlation of [CCH]/[CH₃OH] with L/M indicates that in the central region of massive clumps this ratio serves as an evolutionary indicator. The ratio [CH₃CN]/[CH₃OH] has a correlation with L/M as well: except the hot massive core G31, the other sources show a clear monotonic increase in [CH₃CN]/[CH₃OH] with increasing L/M. This is consistent with the formation timescale of the two species: CH₃OH forms at early stages (10 K) on grain surfaces, primarily via successive hydrogenation of CO (e.g., Watanabe et al. 2003), while the formation of CH₃CN appears late through radiative association between HCN and CH₃⁺ on the grain surface and between CH₃ and CN in the gas phase (Nomura & Millar 2004).

### 4.4.2. Sulfur-bearing molecules

Time-dependent theoretical chemical models of sulfur-bearing species find that the evaporation of H₂S and subsequent chemical reactions in hot cores lead to an increasing abundance ratio of SO₂/SO (Wakelam et al. 2011; Esplugues et al. 2014). At a late stage of hot core evolution, H₂CS and CS are also drastically enhanced; together with SO₂, they become the most abundant sulfur-bearing species. In shocked regions similar trends along evolution are predicted (Wakelam et al. 2005; Esplugues et al. 2014), but observations seem to suggest opposite variations of X(CS)/X(SO) with source evolutionary stages (Li et al. 2015; Gerner et al. 2014; Herpin et al. 2009). Assuming the same ³⁴S/S isotopic ratio for our target clumps, we find a similar result, which is reflected by the feature in Fig. 19 (left panel, bottom row). Similar results of X(H₂CS)/X(SO) and X(SO₂)/X(SO) with clump L/M are derived in Fig. 19 (middle and right panel, bottom row); there are no prominent correlations found. We note that the temporal evolution of abundance of sulfur-bearing species is highly dependent on the source physical structure (Wakelam et al. 2005; Wakelam et al. 2011), which require

tailored chemical modeling for individual source to disentangle this impact and the evolutionary time.

### 4.5. Density structure: relation with cloud structure and implications for different dense gas conversions

In this section we discuss how the radial density profiles of massive clumps (see Sect. 4.2) may be linked to statistics of parental cloud structures on large scales (≥10 pc). We also elaborate on the physical implications of the radial density profiles comparisons between bulk gas and dense gas structures, $\rho_{\text{bulk}}$ and $\rho_{\text{dense}}$.

Gravoturbulence simulations of molecular clouds (e.g., Kritsuk et al. 2011; Lee et al. 2015) showed that the cloud volume density probability distribution functions ($\rho$-PDF, $p_s$) have lognormal and power-law forms in the low- and high-density regimes, respectively (Fig. 1). The lognormal form is attributed to the primordial and maintained supersonic turbulence, while the power-law tail is usually regarded as a sign of gravitational collapse. Thus, the study of PDF statistics is a useful method to understand the physical mechanisms at play in molecular clouds. Gravitational contraction converts a fraction of low-density gas to high-density structures over a timescale comparable to the free-fall timescale ($t_{\text{ff}}$) of the mean density. With this process, the slope of the power-law tail $s$ changes from $s \approx -3$ to $s \approx -1.5-1$ (e.g., Girichidis et al. 2014; Guszejnov et al. 2018) in a mean free-fall time. When the power-law tail of the $\rho$-PDF is dominated by a single gravitationally bound massive clump with a $\rho(r) \propto r^q$ density profile, there is a relation $s=3/q$ (Federrath & Klessen 2013). In other words, $q$ is expected to evolve from $\approx -1$ to $-3$ to $-2$ over a mean $t_{\text{ff}}$. Observational works usually constrain the column density probability distribution function (N-PDF, $p_{\eta}$) rather than $\rho$-PDF. In this case the power-law slope of the N-PDF ($s_N$) can be related to the radial density profile of the molecular clump by $\frac{2}{1+q}$. Column density mapping toward Galactic massive star-forming complexes yield power-law slopes for N-PDF ranging between $\sim -4$ to $-2$ (Lin et al. 2017). The evolutionary trend of steepening density profiles (the bulk gas density profile) of the massive clumps resolved in this work is in general consistent with the prediction from cloud-scale statistics.

The cloud-scale (10 pc) dense gas fraction can be relatively well described by statistics of $\rho/N$-PDFs, which do not factor in features of spatial distribution. Within the $\sim 1$ pc scale clump structure, the spatial configuration of dense gas can be critical in determining the properties of final star clusters. We discussed the gas radial density profiles probed by the dust emission and dense gas tracer (CH₃OH) separately in Sect. 4.2. Here we note again that dust emission traces averaged (mass-weighted) bulk gas structures ($\bar{\rho} \sim 10^{4.5}$ cm⁻³), while CH₃OH emission probes dense gas regimes ($\rho \gtrsim 10^{6.5}$ cm⁻³, up to several $10^8$ cm⁻³) within the massive clumps. We note that the gas density difference seen by dust and higher angular resolution CH₃OH emission can also be partially related to the coarser resolution of the single-dish continuum images, which tend to smear out the gas overdensities at smaller spatial scales. However, the major cause of the difference is that CH₃OH transitions selectively trace the higher density gas component, while the dust emission does not critically depend on the volume density of the gas (in the regimes we are interested in), providing a measure of bulk gas properties. By relying on two density probes for different gas density regimes, we go beyond the self-similar solution of gravitational collapse by characterizing the density configuration of massive clumps with two layers, $\rho_{\text{bulk}}$ and $\rho_{\text{dense}}$, each of which is described by a power-law form. The $\rho_{\text{dense}}$ traced by CH₃OH from RADEX
<figure><img src="image_24.png" /><figcaption>Fig. 22. Comparison between gas density profiles derived by modeling continuum and CH3OH line emission. The ratio of densities derived by CH3OH LIME modeling ($\rho_{dense}$) to continuum results ($\rho_{bulk}$) are shown as pluses (following right y-axis). The density estimated by SMA 1.2 mm continuum observations representing the central core average density is shown as a vertical orange line in each plot (Table 4).</figcaption></figure>

modeling (Sect. 3.5) is verified by full radiative transfer non-LTE models described in Appendix F, by introducing a volume filling factor $ff_{dens}$ (Table 7) to spherical models. The low volume filling factor indicates that the dense gas is highly clumpy and may be organized into filament- or sheet-like structures, which is a natural outcome of supersonic turbulence and gravitational collapse (Smith et al. 2020). In general there seems to be a higher density contrast (decreasing $ff_{dens}$) with increasing $L/M$ (listed sequentially in Table 7, last column). While $ff_{dens}$ provides an upper limit of the volume filling factor of high-density gas, we make some inferences here for the dense gas mass fraction based on intuitive assumptions and references to numerical simulations. The aim is to establish the ratio of the bulk gas to dense gas density profiles, $\rho_{bulk}/\rho_{dense}$, as a measure of the dense gas fraction, and the radial profile of this ratio as an indicator of different star formation modes (or temporal variations of the SFE) for the massive clumps.

It is a well-known fact that there is a strong correlation between SFR and dense gas. Based on a gravo-turbulent fragmentation scenario, Padoan et al. (2012, 2017) derived an empirical relation between SFR per free-fall time $\epsilon_{ff}$ (the fraction of a cloud’s mass that is converted to stellar mass over a free-fall timescale, Krumholz & McKee 2005) and $\alpha_{vir}$,

$$ \epsilon_{ff} = \dot{M}/(M/t_{ff}) \propto \exp(-\alpha_{vir}^{1/2}). \qquad (8) $$

In a hierarchical description of cloud structure where gas density increases with decreasing scale, the conservation of $\dot{M}$ (or equivalently SFR) translates into

$$ \epsilon_{ff,1} M_{1,tot} / t_{ff,\bar{\phi}_1} = \epsilon_{ff,2} M_{2,tot} / t_{ff,\bar{\phi}_2}, \qquad (9) $$

where the number indices of 1 and 2 denote two adjacent levels in a hierarchy. We assume $\rho_{bulk}$ represents the gas density of a lower level structure in the hierarchy from which a higher level structure of gas density $\rho_{dense}$ originates from. We further assume that only the dense gas $\rho_{dense}$ participates in the star formation process. Then with a scale-invariant $\epsilon_{ff}$ (or equivalently scale-invariant $\alpha_{vir}$) Eq. (9) translates to $M_{dense} = M_{tot}(t_{ff,\rho_{dense}}/t_{ff,\rho_{bulk}})$.

Substituting the free-fall timescale $t_{ff} \propto \rho^{-1/2}$ into the equation, we obtain

$$ M_{dense} \sim M_{tot} \left( \frac{\rho_{bulk}}{\rho_{dense}} \right)^{1/2}. \qquad (10) $$

Naively, this relation means that at a certain scale, the closer the derived dense gas density is to the bulk gas density, the more mass fraction of the clump is occupied by the dense gas. We note that even without the assumption of scale-invariant $\epsilon_{ff}$, considering only that the two set of densities are each the mass-weighted value in their respective density regimes, the mass fraction of dense gas to total clump gas is also largely determined by $\rho_{bulk}/\rho_{dense}$ when the ratio between density of bulk gas and gas not seen by CH3OH is similar among sources and does not vary with scale.

Figure 22 shows the density profiles $\rho_{dense}$ and $\rho_{bulk}$. This figure also shows the ratios between these two profiles, which provide a measure of dense gas mass fraction (DGMF): $M_{dense}/M_{tot}$. We found that the DGMF in G19, G08b, and G10 decreases with radius, while it increases with radius in G28, G08a, G13, and G31. For G31, the optical depths of the CH3OH lines are high such that the line ratios are less effective densitometers in the density regimes it associates with than the other sources. We thus omit discussion on this source, while noting that the clump indeed holds the highest density values for dense gas seen by CH3OH among all target sources. Of the other six sources, G19 shows a decreasing DGMF from 25% to 19%, G08b from 22% to 18%, and G10 from 24% to 15%. For G28, the DGMF increases from 15% to ~17% at 0.1–0.2 pc. In G08a and G13, the DGMF is ~19% at their centers and achieve ~22% at 0.2–0.3 pc radii.

Given these radial variations of DGMF, it might be indicative that of these six sources, G19 and G10 may be most efficiently converting gas into the dense gas regime, having a focused DGMF toward clump center. G08b exhibits a similarly high DGMF that changes slightly with scale, while early-phase clump G08a and the least massive G13 have DGMF peaking around the intermediate scales relative to the center. Another early-stage
source, G28, has the smallest DGMF overall and also shows a relatively invariant DGMF over scales. This trend is only partially reflected in the density profile slope of $\rho_{\text{bulk}}(r)$, which is generally regarded as a measure of dense gas concentration that directly inflates the SFR (e.g., Parmentier 2019). It is noteworthy that radial variations of the DGMF for individual sources are in general compatible with the radial dependence of the virial parameter shown in Sect. 4.3 (Fig. 21). These two measures are both indicative of the capability of the source in converting gas into the denser regime at different positions relative to the bottom of the clump’s gravitational potential. It is clear that the three classes of increasing, flattened, and decreasing DGMF as a function of radius, for a particular source, correspond respectively to a shallower, similar, and steeper dense gas profile ($q_{\text{dense}}$) compared to that of its bulk gas ($q_{\text{bulk}}$) on the basis of Eqs. (2), (4), and 10. In the non-homologous spherical collapse framework, the physical condition of an outer shell can be regarded as an earlier stage preceding the inner shell (e.g., Vázquez-Semadeni et al. 2019), then a radially decreasing (increasing) DGMF corresponds to an accelerating (retarding) dense gas conversion and temporal increase (decrease) of SFE for a particular source.

The localized behavior of dense gas conversion may be causing the chaotic and scattered SFR versus gas density relation, and this behavior goes beyond the self-similar solutions for collapsing clumps. The non-self-similar behavior is suggested to be relevant even for the spherical collapse, as shown theoretically by Coughlin (2017), dependent on the initial conditions. The spatial variations of the DGMF may also reflect two competing gravitational collapses within these massive clumps: collapse toward the global potential center and collapse of dense regions into ambient filaments surrounding the center. The latter process is mostly induced by turbulence (Girichidis et al. 2011). The dominance of one or another process may depend on the initial density profile and turbulence driving modes (Girichidis et al. 2011; Lomax et al. 2015).

We can also establish $\rho_{\text{bulk}}/\rho_{\text{dense}}$ as a measure of DGMF from another perspective, with a weaker assumption. In the framework of turbulent convergent flows, the gas density enhancement after a compressive shock has an inverse length scale relation following

$$ \rho_{\text{post}}/\rho_0 \propto L_0/l_{\text{post}} \sim M_s^2, \quad (11) $$

where $L_0$ and $l_{\text{post}}$ denote the length scales of the preshock and postshock gas, and $\rho_0$ and $\rho_{\text{post}}$ the gas densities, respectively; $M_s$ denotes the sonic Mach number, $M_s = \frac{\sigma_{\text{rms}}}{c_s}$. This shock jump condition is at the origin of the lognormal distribution of PDF. As Eq. (11) links the gas length scale to density with an inverse relation, then the mass ratio of postshock dense gas versus the preshock gas may also be represented by a power-law form following $(\rho_0/\rho_{\text{post}})^{-s}$ with slope $s$ dependent on the assumed geometry. With a cylindrical geometry describing infinite filaments, the gas mass is $\propto \rho l^2$, where $l$ denotes thickness (radius), with slope $s=1$.

If the enhancement of gas density ($\rho_{\text{dense}}$) as traced by CH$_3$OH is regarded as the result of compressive turbulence on the preshock gas that has a density represented by $\rho_{\text{bulk}}$, the ratio $(\rho_0/\rho_{\text{post}})$ still holds as a DGMF measure. Compared to Eq. (10) ($s=1/2$), the different scaling of $s=1$ does not change the general trend of the radial change of this measure, but only increases the contrast between the outer layer and inner region.

In any case, the dense gas probed by CH$_3$OH does have contributions from shock entrainment, as suggested by locally

increased $M_s$ and the fact that CH$_3$OH is likely enhanced in shocked regions. On the other hand, we have already seen in Sect. 4.3 that gravitational collapse has altered the general scaling relation of $v_T$ (hence $M_s$), which makes the slope $s$ vary. A more stringent comparison of radial change of dense gas mass fraction to indicate the dense gas conversion efficiency would benefit from properly separating the part of dense gas associated with transient gas substructures with virialized cores and coherent flows, which is achievable with finer resolution observations (both spectral and spatial).

## 5. Conclusions

The gas thermal properties are critical to the star formation process. For massive star formation, the density structure of the nursery gas clumps may not be well predicted by previous hydrostatic equilibrium core models, which are generally representative of low-mass cores, due to a higher level of turbulence and gravitational collapse set at larger scales. We conducted a SMA and APEX line survey toward a sample of eight massive star-forming clumps in order to understand the evolution of temperature and density structures. The major findings are as follows:

1. Transitions of multiple molecular species (CH$_3$CN, CH$_3$CCH, H$_2$CS) of distinct critical densities, together with dust emission, provide a reasonably good sampling of different gas temperature regimes ($\ge 200-20$ K) over the full clump scale. There is not a single power-law relation dependent only on source luminosity that could describe all the radial temperature variations in our sample. The elevated temperature of a less luminous source may be related to the intermittency of accretion and shock-related activities.

2. CH$_3$OH line series are good density probes for massive clumps, selectively tracing density regimes of $\ge 10^6$ cm$^{-3}$. Systematic steepening of density profiles along clump evolution, indicated by $L/M$, is revealed in the sample via continuum (bulk gas) and CH$_3$OH line (dense gas) modeling. The density slopes change from $> -1$ to $\sim -1.5$. The dense gas component becomes denser with increasing $L/M$, and the volume filling factor of which decreases, corresponding to a higher density contrast along clump evolution.

3. The radial linewidth profile traced by multiple lines displays a scale-dependent relation. Several sources, which have comparable stellar mass and gas content at the scale probed, have a central decreasing trend of $v_T \propto r^{-0.4}$ and $v_T \propto r^{0.2}$ in the outer envelope, which may be related to conversion of gravitational energy to turbulence. On larger scales (0.2–0.5 pc) all clumps are close to virial state ($\alpha_{\text{vir}}/\alpha_{\text{cri}} \sim 1$), as indicated by CH$_3$CCH, H$^{13}$CO+ lines. Small scale ($\le 0.1$ pc) virial parameters traced by CH$_3$CN and H$_2$CS can exceed equilibrium values for some of the sources, the origin of which is hindered by our limited resolution. Overall, we observe that the clump gas evolves from a super-virial to sub-virial state with increasing $L/M$.

4. The abundances of CH$_3$OH, CH$_3$CN, and H$_2$CS show better correlation with source central temperature than with luminosity. Abundance ratios of [CCH]/[CH$_3$OH] and [CH$_3$CN]/[CH$_3$OH] are in good correlation with clump $L/M$, and can be used as indicators of evolutionary stages of massive star-forming clumps.

5. The evolutionary trend of clump density profiles is compatible with cloud-scale diagnosis that frequently reveal the time-varying power-law tail of PDFs. In a hierarchical view, the radial variation in the dense gas mass fraction, which can
be approximated by the density ratio representing the aver-
aged bulk gas (dust continuum, no spatial filtering) to that
probed by a high-density tracer (CH3OH), may be indica-
tive of the efficiency of the source in dense gas conversion
and dense gas focusing to the global gravitational center. The
bulk gas density profile is a less distinct measure because
self-similar approximation does not adequately describe the
clumpy gas structures within massive clumps.

In this work, we attempted to understand the role of self-gravity
in shaping the gas structure as massive clumps evolve over
time. Our results are based on ≤0.1 pc observations, although a
dense gas filling factor (probed by CH3OH) indicates, on smaller
scales, a clumpy gas environment of higher density contrast
along the evolutionary track. Detailed fragmentation properties,
cores and filaments, and their kinematics remain unresolved.
In any case, deduction of line-of-sight cloud geometry is diffi-
cult. With the unknown cloud thickness, the mass scale of the
high-density gas regime remains highly uncertain.

Existing observations have broadly revealed transonic to
subsonic linewidths associated with localized substructures in
massive star-forming regions, the length scale of which is well
above sonic scale, indicating an efficient turbulence dissipata-
tion process. More robust velocity linewidth profiles require
better resolved observations to disentangle multiple velocity
components, and confusion from stellar feedback (outflows). In
addition, while the scale-dependent linewidth versus radius rela-
tion is indirect evidence (based on the assumption that it mirrors
radial infall velocity) of the changing role of gravitational col-
lapse in dominating gas dynamics, more direct evidence would
be measurement of infall velocities, or accretion rates as a func-
tion of radius. The latter is indispensable to distinguish different
collapse models (e.g., Padoan et al. 2020). This calls for obser-
vations of multiple “infall tracers” sensitive to different gas
density regimes, along with a proper description of the physical
structure.

Acknowledgements. We are thankful to the referee whose comments and suggestions helped us to greatly increase the clarity of the paper. This work was partly funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Collaborative Research Council 956, sub-project A6. Y.L. is a member of the International Max-Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. H.B.L. is supported by the Ministry of Science and Technology (MoST) of Taiwan (Grant Nos. 108-2112-M-001-002-MY3 and 110-2112-M-001-069-). The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.

References

Adams, F. C. 1991, ApJ, 382, 544

Adams, F. C., & Shu, F. H. 1985, ApJ, 296, 655

Bachiller, R., & Gutiérrez, M. P. 1997, ApJ, 487, L93

Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6

Ballesteros-Paredes, J., Vázquez-Semadeni, E., Gazol, A., et al. 2011, MNRAS, 416, 1436

Ballesteros-Paredes, J., Vázquez-Semadeni, E., Palau, A., & Klessen, R. S. 2018, MNRAS, 479, 2112

Bate, M. R. 2009, MNRAS, 392, 1363

Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250

Bell, T. A., Cernicharo, J., Viti, S., et al. 2014, A&A, 564, A114

Belloche, A., Meshcheryakov, A. A., Garrod, R. T., et al. 2017, A&A, 601, A49

Beltrán, M. T., Cesaroni, R., Neri, R., et al. 2004, ApJ, 601, L187

Beltrán, M. T., Cesaroni, R., Rivilla, V. M., et al. 2018, A&A, 615, A141

Bergin, E. A., Goldsmith, P. F., Snell, R. L., & Ungerechts, H. 1994, ApJ, 431, 674

Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140

Beuther, H., Schilke, P., Menten, K. M., et al. 2002, ApJ, 566, 945

Beuther, H., Leurini, S., Schilke, P., et al. 2007, A&A, 466, 1065

Beuther, H., Semenov, D., Henning, T., & Linz, H. 2008, ApJ, 675, L33

Beuther, H., Henning, T., Linz, H., et al. 2010, A&A, 518, L78

Beuther, H., Ragan, S. E., Johnston, K., et al. 2015, A&A, 584, A67

Binder, B. A., & Povich, M. S. 2018, ApJ, 864, 136

Bisschop, S. E., Jørgensen, J. K., van Dishoeck, E. F., & de Wachter, E. B. M. 2007, A&A, 465, 913

Blake, G. A., van Dishoeck, E. F., Jansen, D. J., Groesbeck, T. D., & Mundy, L. G. 1994, ApJ, 428, 680

Bøgelund, E. G., Barr, A. G., Taquet, V., et al. 2019, A&A, 628, A2

Bracco, A., Palmeirim, P., André, P., et al. 2017, A&A, 604, A52

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2021, astropy/photutils: 1.2.0

Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25

Calcutt, H., Willis, E. R., Jørgensen, J. K., et al. 2019, A&A, 631, A137

Camacho, V., Vázquez-Semadeni, E., Ballesteros-Paredes, J., et al. 2016, ApJ, 833, 113

Camacho, V., Vázquez-Semadeni, E., Palau, A., Busquet, G., & Zamora-Avilés, M. 2020, ApJ, 903, 46

Carey, S. J., Feldman, P. A., Redman, R. O., et al. 2000, ApJ, 543, L157

Carey, S. J., Noriega-Crespo, A., Mizuno, D. R., et al. 2009, PASP, 121, 76

Caselli, P., & Myers, P. C. 1995, ApJ, 446, 665

Caswell, J. L. 1998, MNRAS, 297, 215

Cesaroni, R., Olmi, L., Walmsley, C. M., Churchwell, E., & Hofner, P. 1994, ApJ, 435, L137

Cesaroni, R., Beltrán, M. T., Zhang, Q., Beuther, H., & Fallscheer, C. 2011, A&A, 533, A73

Cesaroni, R., Beuther, H., Ahmadi, A., et al. 2019, A&A, 627, A68

Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Churchwell, E., Walmsley, C. M., & Cesaroni, R. 1990, A&AS, 83, 119

Contreras, Y., Schuller, F., Urquhart, J. S., et al. 2013, A&A, 549, A45

Coughlin, E. R. 2017, ApJ, 835, 40

Csengeri, T., Bontemps, S., Schneider, N., et al. 2011, ApJ, 740, L5

Csengeri, T., Urquhart, J. S., Schuller, F., et al. 2014, A&A, 565, A75

Csengeri, T., Weiss, A., Wyrowski, F., et al. 2016, A&A, 585, A104

Cummins, S. E., Green, S., Thaddeus, P., & Linke, R. A. 1983, ApJ, 266, 331

Curry, C. L., & McKee, C. F. 2000, ApJ, 528, 734

Cyganowski, C. J., Whitney, B. A., Holden, E., et al. 2008, AJ, 136, 2391

Dempsey, J. T., Friberg, P., Jenness, T., et al. 2013, MNRAS, 430, 2534

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015]

Esplugues, G. B., Viti, S., Goicoechea, J. R., & Cernicharo, J. 2014, A&A, 567, A95

Estalella, R. 2017, PASP, 129, 025003

Falgarone, E., & Puget, J. L. 1985, A&A, 142, 157

Fallscheer, C., Beuther, H., Zhang, Q., Keto, E., & Sridharan, T. K. 2009, A&A, 504, 127

Fayolle, E. C., Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Bisschop, S. E. 2015, A&A, 576, A45

Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621

Flower, D. R., Pineau des Forêts, G., & Rabli, D. 2010, MNRAS, 409, 29

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303

Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927

Garrod, R. T., Belloche, A., Müller, H. S. P., & Menten, K. M. 2017, A&A, 601, A48

Gerner, T., Beuther, H., Semenov, D., et al. 2014, A&A, 563, A97

Giannetti, A., Leurini, S., Wyrowski, F., et al. 2017, A&A, 603, A33

Gieser, C., Beuther, H., Semenov, D., et al. 2021, A&A, 648, A66

Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011, MNRAS, 413, 2741

Girichidis, P., Konstandin, L., Whitworth, A. P., & Klessen, R. S. 2014, ApJ, 781, 91

Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 421, 9

Goldsmith, P. F. 2001, ApJ, 557, 736

Gómez, L., Luis, L., Hernández-Curiel, I., et al. 2010, ApJS, 191, 207

Gómez, G. C., Vázquez-Semadeni, E., & Palau, A. 2021, MNRAS, 502, 4963

Gómez-Ruiz, A. I., Kurtz, S. E., Araya, E. D., Hofner, P., & Loinard, L. 2016, ApJS, 222, 18

Graninger, D. M., Wilkins, O. H., & Öberg, K. I. 2016, ApJ, 819, 140

Grave, J. M. C., & Kumar, M. S. N. 2009, A&A, 498, 147

Green, S. 1986, ApJ, 309, 331

Green, J. A., & McClure-Griffiths, N. M. 2011, MNRAS, 417, 2500

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3
Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13

Guszejnov, D., Krumholz, M. R., & Hopkins, P. F. 2016, MNRAS, 458, 673

Guszejnov, D., Hopkins, P. F., & Ma, X. 2017, MNRAS, 472, 2107

Guszejnov, D., Hopkins, P. F., & Grudić, M. Y. 2018, MNRAS, 477, 5139

Hacar, A., Kainulainen, J., Tafalla, M., Beuther, H., & Alves, J. 2016, A&A, 587, A97

Harju, J., Lehtinen, K., Booth, R. S., & Zinchenko, I. 1998, A&AS, 132, 211

Hatchell, J., & van der Tak, F. F. S. 2003, A&A, 409, 589

Helmich, F. P., & van Dishoeck, E. F. 1997, A&AS, 124, 205

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020, ApJ, 904, 194

Henriksen, R., Andre, P., & Bontemps, S. 1997, A&A, 323, 549

Herpin, F., Marseille, M., Wakelam, V., Bontemps, S., & Lis, D. C. 2009, A&A, 504, 853

Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092

Hildebrand, R. H. 1983, QJRAS, 24, 267

Ho, P. T. P., & Haschick, A. D. 1986, ApJ, 304, 501

Hofner, P., & Churchwell, E. 1996, A&AS, 120, 283

Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823

Humire, P. K., Thiel, V., Henkel, C., et al. 2020, A&A, 642, A222

Izquierdo, A. F., Galván-Madrid, R., Maud, L. T., et al. 2018, MNRAS, 478, 2505

Izquierdo, A. F., Smith, R. J., Glover, S. C. O., et al. 2021, MNRAS, 500, 5268

Jacobsen, S. K., Jørgensen, J. K., van der Wiel, M. H. D., et al. 2018, A&A, 612, A72

Jappsen, A. K., Klessen, R. S., Larson, R. B., Li, Y., & Mac Low, M. M. 2005, A&A, 435, 611

Jiang, X.-J., Liu, H. B., Zhang, Q., et al. 2015, ApJ, 808, 114

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007, ApJ, 659, 479

Juvela, M., He, J., Pattle, K., et al. 2018, A&A, 612, A71

Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35

Kainulainen, J., Beuther, H., Banerjee, R., Federrath, C., & Henning, T. 2011, A&A, 530, A64

Kauffmann, J., Pillai, T., & Goldsmith, P. F. 2013, ApJ, 779, 185

Keto, E. R. 1990, ApJ, 355, 190

Keto, E. R., Ho, P. T. P., & Reid, M. J. 1987, ApJ, 323, L171

Keto, E. R., Ho, P. T. P., & Haschick, A. D. 1988, ApJ, 324, 920

Klessen, R. S. 2000, ApJ, 535, 869

Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17

Klessen, R. S., Spaans, M., & Jappsen, A.-K. 2007, MNRAS, 374, L29

Koda, J., Sawada, T., Wright, M. C. H., et al. 2011, ApJS, 193, 19

Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20

Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545

Kroupa, P., Weidner, C., Pflamm-Altenburg, J., et al. 2013, The Stellar and Sub- Stellar Initial Mass Function of Simple and Composite Populations, eds. T. D. Oswalt & G. Gilmore, 5, 115

Krumholz, M. R. 2014, MNRAS, 437, 1662

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Krumholz, M. R., Klein, R. I. & McKee, C. F. 2015, ApJ, 656, 959

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71

Kuiper, T. B. H., Kuiper, E. N. R., Dickinson, D. F., Turner, B. E., & Zuckerman, B. 1984, ApJ, 276, 211

Kurono, Y., Morita, K.-I., & Kamazaki, T. 2009, PASJ, 61, 873

Larson, R. B. 1969, MNRAS, 145, 271

Larson, R. B. 1981, MNRAS, 194, 809

Lee, Y.-N., & Hennebelle, P. 2018, A&A, 611, A88

Lee, Y.-N., & Hennebelle, P. 2019, A&A, 622, A125

Lee, E. J., Chang, P., & Murray, N. 2015, ApJ, 800, 49

Leurini, S., Schilke, P., Menten, K. M., et al. 2004, A&A, 422, 573

Leurini, S., Schilke, P., Wyrowski, F., & Menten, K. M. 2007, A&A, 466, 215

Leurini, S., Codella, C., López-Sepulcre, A., et al. 2014, A&A, 570, A49

Leurini, S., Menten, K. M., & Walmsley, C. M. 2016, A&A, 592, A31

Li, J., Wang, J., Zhu, Q., Zhang, J., & Li, D. 2015, ApJ, 802, 40

Li, J. I.-H., Liu, H. B., Hasegawa, Y., & Hirano, N. 2017, ApJ, 840, 72

Li, S., Zhang, Q., Pillai, T., et al. 2019, ApJ, 886, 130

Lin, Y., Liu, H. B., Li, D., et al. 2016, ApJ, 828, 32

Lin, Y., Liu, H. B., Dale, J. E., et al. 2017, ApJ, 840, 22

Lin, Y., Csengeri, T., Wyrowski, F., et al. 2019, A&A, 631, A72

Liu, H. B. 2017, A&A, 597, A70

Liu, H. B. 2019, ApJ, 877, L22

Liu, H. B., Ho, P. T. P., Zhang, Q., et al. 2010, ApJ, 722, 262

Liu, H. B., Zhang, Q., & Ho, P. T. P. 2011, ApJ, 729, 100

Liu, H. B., Quintana-Lacaci, G., Wang, K., et al. 2012, ApJ, 745, 61

Liu, H. B., Galván-Madrid, R., Jiménez-Serra, I., et al. 2015, ApJ, 804, 37

Liu, H. B., Chen, H.-R. V., Román-Zúñiga, C. G., et al. 2019, ApJ, 871, 185

Lodders, K. 2003, ApJ, 591, 1220

Lomax, O., Whitworth, A. P., & Hubber, D. A. 2015, MNRAS, 449, 662

Longmore, S. N., Pillai, T., Keto, E., Zhang, Q., & Qiu, K. 2011, ApJ, 726, 97

Lowe, I., Mason, B., Bhandarkar, T., et al. 2021, ArXiv e-prints, [arXiv:2105.13432]

Lynch, S. M. 2007, Introduction to Applied Bayesian Statistics and Estimation for Social Scientists (New York, NY: Springer Science & Business Media, LLC.)

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

McKee, C. F., & Holliman, John H., I. 1999, ApJ, 522, 313

McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850

McLaughlin, D. E., & Pudritz, R. E. 1997, ApJ, 476, 750

Minh, Y. C., Liu, S. Y., Chen, H. R., & Su, Y. N. 2011, ApJ, 737, L25

Molinari, S., Pezzuto, S., Cesaroni, R., et al. 2008, A&A, 481, 345

Molinari, S., Swinyard, B., Bally, J., et al. 2010, A&A, 518, L100

Molinari, S., Merello, M., Elia, D., et al. 2016, ApJ, 826, L8

Möller, T., Endres, C., & Schilke, P. 2017, A&A, 598, A7

Monsch, K., Pineda, J. E., Liu, H. B., et al. 2018, ApJ, 861, 77

Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41

Mottram, J. C., Beuther, H., Ahmadi, A., et al. 2020, A&A, 636, A118

Mouschovias, T. C., & Morton, S. A. 1991, ApJ, 371, 296

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49

Mueller, K. E., Shirley, Y. L., Evans, Neal J., I., & Jacobson, H. R. 2002, ApJS, 143, 469

Murray, N., & Chang, P. 2015, ApJ, 804, 44

Murray, D. W., Chang, P., Murray, N. W., & Pittman, J. 2017, MNRAS, 465, 1316

Myers, P. C. 2015, ApJ, 806, 226

Nomura, H., & Millar, T. J. 2004, A&A, 414, 409

Öberg, K. I., Fayolle, E. C., Reiter, J. B., & Cyganowski, K. 2014, Faraday Discuss., 168, 81

Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJ, 759, L27

Padoan, P., Haugbølle, T., Nordlund, Å., & Frimann, S. 2017, ApJ, 840, 48

Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82

Palau, A., Estalella, R., Girart, J. M., et al. 2014, ApJ, 785, 42

Pan, L., & Padoan, P. 2009, ApJ, 692, 594

Parmentier, G. 2019, ApJ, 887, 179

Passot, T., & Vázquez-Semadeni, E. 1998, Phys. Rev. E, 58, 4501

Penston, M. V. 1969, MNRAS, 144, 425

Persson, M. V., Harsono, D., Tobin, J. J., et al. 2016, A&A, 590, A33

Pham, D. T., Ghanbarzadeh, A., Koç, E., et al. 2006, in The Bees Algorithm - A Novel Tool for Complex Optimisation Problems

Pillai, T., Wyrowski, F., Carey, S. J., & Menten, K. M. 2006, A&A, 450, 569

Plume, R., Jaffe, D. T., II, N. J. E., Martin-Pintado, J., & Gomez-Gonzalez, J. 1997, ApJ, 476, 730

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Purcell, C. R., Balasubramanyam, R., Burton, M. G., et al. 2006, MNRAS, 367, 553

Qi, C. 2003, in SFChem 2002: Chemistry as a Diagnostic of Star Formation, eds. C. L. Curry, & M. Fich, 393

Qiu, K., Zhang, Q., Beuther, H., & Yang, J. 2007, ApJ, 654, 361

Rabli, D., & Flower, D. R. 2010, MNRAS, 406, 95

Ragan, S., Henning, T., Krause, O., et al. 2012, A&A, 547, A49

Ragan, S. E., Henning, T., Tackenberg, J., et al. 2014, A&A, 568, A73

Rivilla, V. M., Beltrán, M. T., Cesaroni, R., et al. 2017, A&A, 598, A59

Robertson, B., & Goldreich, P. 2012, ApJ, 750, L31

Rodríguez-Garza, C. B., Kurtz, S. E., Gómez-Ruiz, A. I., et al. 2017, ApJS, 233, 4

Rolffs, R., Schilke, P., Wyrowski, F., et al. 2011, A&A, 527, A68

Roy, A., André, P., Palmeirim, P., et al. 2014, A&A, 562, A138

Sanhueza, P., Jackson, J. M., Zhang, Q., et al. 2017, ApJ, 841, 97

Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102

Sanna, A., Reid, M. J., Menten, K. M., et al. 2014, ApJ, 781, 108

Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in ASP Conf. Ser., 77, Astronomical Data Analysis Software and Systems IV, eds. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 433

Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415

Shirley, Y. L. 2015, PASP, 127, 299
Shu, F. H. 1977, ApJ, 214, 488

Siringo, G., Kreysa, E., Kovács, A., et al. 2009, A&A, 497, 945

Smith, R. J., Longmore, S., & Bonnell, I. 2009, MNRAS, 400, 1775

Smith, R. J., Treß, R. G., Sormani, M. C., et al. 2020, MNRAS, 492, 1594

Sollins, P. K., Zhang, Q., Keto, E., & Ho, P. T. P. 2005, ApJ, 624, L49

Spaans, M., & Silk, J. 2000, ApJ, 538, 115

Sridharan, T. K., Beuther, H., Saito, M., Wyrowski, F., & Schilke, P. 2005, ApJ, 634, L57

Sutton, E. C., Blake, G. A., Genzel, R., Masson, C. R., & Phillips, T. G. 1986, ApJ, 311, 921

Traficante, A., Lee, Y. N., Hennebelle, P., et al. 2018, A&A, 619, L7

Troscompt, N., Faure, A., Maret, S., et al. 2009, A&A, 506, 1243

Urban, A., Martel, H., & Evans, N. J. II 2010, ApJ, 710, 1343

Urquhart, J. S., König, C., Giannetti, A., et al. 2018, MNRAS, 473, 1059

van der Tak, F. F. S., van Dishoeck, E. F., Evans, Neal J., I., & Blake, G. A. 2000, ApJ, 537, 283

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van Dishoeck, E. F., Blake, G. A., Jansen, D. J., & Groesbeck, T. D. 1995, ApJ, 447, 760

van 't Hoff, M. L. R., van Dishoeck, E. F., Jørgensen, J. K., & Calcutt, H. 2020, A&A, 633, A7

Varricatt, W. P., Davis, C. J., Ramsay, S., & Todd, S. P. 2010, MNRAS, 404, 661

Varricatt, W. P., Wouterloot, J. G. A., Ramsay, S. K., & Davis, C. J. 2018, MNRAS, 480, 4231

Vasyunina, T., Vasyunin, A. I., Herbst, E., et al. 2014, ApJ, 780, 85

Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés. M. 2019, MNRAS, 490, 3061

Wakelam, V., Ceccarelli, C., Castets, A., et al. 2005, A&A, 437, 149

Wakelam, V., Hersant, F., & Herpin, F. 2011, A&A, 529, A112

Walch, S., Burkert, A., Whitworth, A., Naab, T., & Gritschneder, M. 2009, MNRAS, 400, 13

Wang, Y., Zhang, Q., Pillai, T., Wyrowski, F., & Wu, Y. 2008, ApJ, 672, L33

Ward-Thompson, D., Motte, F., & Andre, P. 1999, MNRAS, 305, 143

Watanabe, N., Shiraki, T., & Kouchi, A. 2003, ApJ, 588, L121

Whitworth, A. P., Bhattal, A. S., Francis, N., & Watkins, S. J. 1996, MNRAS, 283, 1061

Wienen, M., Wyrowski, F., Menten, K. M., et al. 2015, A&A, 579, A91

Williams, J. P., Blitz, L., & McKee, C. F. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell, 97

Williams, S. J., Fuller, G. A., & Sridharan, T. K. 2005, A&A, 434, 257

Wilson, T. L., Rohlfs, K., & Hüttemeister, S. 2013, Tools of Radio Astronomy

Wood, D. O. S., & Churchwell, E. 1989, ApJS, 69, 831

Wyrowski, F., Güsten, R., Menten, K. M., Wiesemeyer, H., & Klein, B. 2012, A&A, 542, L15

Xu, S., & Lazarian, A. 2020, ApJ, 890, 157

Zapata, L. A., Rodríguez, L. F., Ho, P. T. P., Beuther, H., & Zhang, Q. 2006, AJ, 131, 939

Zhang, Q., Wang, Y., Pillai, T., & Rathborne, J. 2009, ApJ, 696, 268
# Appendix A: Target sources

## Appendix A.1: G19.882-00.534

G19.882-00.534 (IRAS 18264–1152) is classified as an extended green object (EGO, Cyganowski et al. 2008), located at 3.7 kpc with luminosity of >$10^4 L_\odot$. The most prominent features of this source are the high-velocity outflow (Qiu et al. 2007, Leurini et al. 2014) and the high level of maser activity (among a sample of 56 high-mass star-forming regions, Rodríguez-Garza et al. 2017). The outflow is oriented in an east–west direction, showing enhanced H$_2$ near-infrared emission as well (Varricatt et al. 2010). Sensitive 7 mm and 1.3 cm observations resolved the source into a triple system, consisting of two optically thin HII or dust emitting sources, and a thermal jet or a partially optically thick HII region (Zapata et al. 2006).

## Appendix A.2: G08.684-00.366 and G08.671-00.356

These two massive star-forming clumps are part of the IRAS 18032-2137 star-forming cloud located at 4.8 kpc (Purcell et al. 2006). Multiple water, methanol, and hydroxide masers are detected toward both sources (Hofner & Churchwell 1996, Caswell 1998, Gómez et al. 2010). G08.671-00.356 (~9000 $L_\odot$) is a UCHII region (Wood & Churchwell 1989), while G08.684-00.366 (~3000 $L_\odot$) is a less evolved relatively infrared weak source, lying 1' offset in the northeastern direction. Strong SiO emission toward G08.684-00.366 indicates that star formation already takes place in this source (Harju et al. 1998). Previous SMA observations resolve the source into three dense cores accompanied by extended outflow components traced by CO (2-1) emission; the core masses are in the range of ~2-15 $M_\odot$, a small fraction of the total clump mass (Longmore et al. 2011).

## Appendix A.3: G10.624-00.380

G10.624-00.380 is an extremely luminous and massive (>$10^5 L_\odot$, ~5000 $M_\odot$) OB cluster-forming clump in the galaxy, located at 4.95 kpc (Sanna et al. 2014). Previous high angular resolution submm and centimeter observations revealed that the clump in the central 0.6 pc is a flattened rotating system (Keto et al. 1987, 1988) where multiple UCHII regions are deeply embedded (Ho & Haschick 1986, Sollins et al. 2005); the clump is fed by the converging flows from ambient filamentary clouds (Liu et al. 2011, 2012). Within the HII region, absorption lines indicate that gas accretion continues despite the ionizing and radiative pressure (Keto 1990). Overall, the clumps seems to be in a global collapsing state with a rotating Toomre-unstable disk-like structure in the center (Liu 2017).

## Appendix A.4: G13.658-00.599

G13.658-00.599 (IRAS 18144-1723) is a molecular clump with a luminosity of >$10^4 L_\odot$ at 3.7 kpc. It is associated with multiple water and methanol masers (Gómez-Ruiz et al. 2016). The H$_2$ emission displays a bowshock feature 18" offset in the west from the center of the IRAS emission, which seems to originate from an extended K-band continuum source associated with the IRAS source (Varricatt et al. 2010). The bowshock feature is surrounded by multiple 44 GHz methanol masers (Gómez-Ruiz et al. 2016), indicating it is likely caused by outflow activity. Deep mid-infrared imaging reveals that the central source hosts two YSOs separated by ~10000 au, at different evolutionary stages; the outflow traced by the CO(3-2) line coincides closely

with the H$_2$ emission, and is likely caused by the younger source in formation (Varricatt et al. 2018).

## Appendix A.5: G31.412+00.307

G31.412+00.307 is a well-known hot massive core that has been extensively studied by both single-dish and interferometry observations (e.g., Cesaroni et al. 1994, Beltrán et al. 2004, 2018, Rivilla et al. 2017), a source showing great chemical richness. It is located at a kinematic distance of 7.9 kpc (Churchwell et al. 1990), and has a luminosity of ≥$10^5 L_\odot$. The central hot core structure is massive and compact (~500 $M_\odot$ of ~8000 au in size, Cesaroni et al. 2011, Beltrán et al. 2004). There is an ultra-compact HII region 5" away from the core (Churchwell et al. 1990). The kinematic features consistently support a rotating core experiencing infall motions (Cesaroni et al. 2011, Wyrowski et al. 2012). Recent ALMA observations of higher angular (~1700 au) resolution suggest that the rotation and infall velocities both increase toward the center, and that the core is composed of a main core of size ~5300 au and a satellite core of much smaller mass (Beltrán et al. 2018). The overall monolithic feature of source G31 makes it an ideal source for understanding the high-mass star formation scenario in light of gas mass origin and evolution.

## Appendix A.6: G18.606-00.074

G18.606-00.074 is a massive infrared-dark clump located in a ~4 pc long filamentary cloud associated with IRAS18223 (Carey et al. 2000, Beuther et al. 2002, Sridharan et al. 2005). The parental molecular cloud is part of an even larger (>50 pc) molecular gas filament (e.g., Kainulainen et al. 2011, Ragan et al. 2014), which is undergoing star formation activity at various evolutionary stages (e.g., Beuther et al. 2002, 2007, 2010, Fallscheer et al. 2009, Ragan et al. 2012). G18.606-00.074 (also known as IRDC18223-3, or as core 11 in Beuther et al. 2015), compared to its ambient core and clump structures, appears to have a larger mass reservoir, showing larger linewidth of $N_2H^+$ (Beuther et al. 2015).

## Appendix A.7: G28.397+00.080

G28.397+00.080 (known as P2 in Wang et al. 2008) is a molecular clump located in a massive (>$10^4 M_\odot$) filamentary infrared dark cloud, G28.34+0.06 (Carey et al. 2000, Pillai et al. 2006, Lin et al. 2017). The central region of the clump is not associated with any near-infrared compact source counterpart, while it appears bright at 24 and 70 µm wavelength (Fig. 4, Wang et al. 2008). It also shows line emission of COMs such as CH$_3$OCH$_3$ and CH$_3$CHO (Zhang et al. 2009, Vasyunina et al. 2014), likely hosting massive protostar(s) which heats the gas up to 45 K, as measured from NH$_3$ observations (Zhang et al. 2009).

# Appendix B: Subtraction of free-free emission from SMA 1.2 mm continuum

We used centimeter radio continuum data collected from the NRAO data archive<fnref n="5" /> for Very Large Array (VLA) for sources G08b, G31 and G10 to estimate the free-free emission contribution to the 1.2 mm flux. Given the time variability of centimeter emission, we selected the most recent data products, whose basic information is listed in Table B.1.

<footnote n="5">⁵ http://www.aoc.nrao.edu/~vlbacald/ArchIndex.shtml</footnote>
Assuming that the dust and free-free emission are both optically thin, and that the dust emission follows a spectral index $\beta = 1.0$, then the relations $S_v \propto v^{2+\beta}$ and $S_v \propto v^{-0.1}$ describe the variation in the two emission components as a function of frequency. The free-free emission at 1.2 mm can then be solved using simultaneous equations of total flux, at two considered frequencies. The centimeter data was smoothed and regridded to match the 1.2 mm data and subtraction was done on a pixel-by-pixel basis. The total flux contributed from free-free emission to the 1.2 mm continuum is also listed in Table B.1.

## Appendix C: RADEX modeling of CH$_3$OH lines and the MCMC procedure

<figure><img src="image_25.png" /><figcaption>Fig. C.1. Posterior distribution of parameters from MCMC RADEX fitting of CH$_3$OH (5-4) line series. The column densities of A/E type are velocity averaged values (log cm$^{-2}$/km s$^{-1}$). The vertical dashed lines in the 1D histograms show the quantiles of 10%, 25%, 50%, 75%, and 90%. The contour levels in the 2D histograms indicate 0.5$\sigma$, 1$\sigma$, 1.5$\sigma$, and 2$\sigma$. The figure shows an example of the fitted parameters of observed lines in one pixel of clump G08b.</figcaption></figure>

For both E-type and A-type CH$_3$OH, we generated RADEX model grids in the column density (N/Δv) range from 10$^{12}$ to 10$^{18}$ cm$^{-2}$ (with 60 logarithmically spaced uniform intervals), in the density range from 10$^4$ to 10$^9$ cm$^{-3}$ (with 100 logarithmically spaced uniform intervals), and in the kinetic temperature range from 10 to 200 K (with 80 uniform intervals). In the RADEX models the linewidth is taken as 1 km/s for all the grids of parameters, while the total column densities (shown in Fig. C.2) are obtained by multiplying the fitted N/Δv with the linewidth obtained from Gaussian fits (Sect. 3.5). The external radiation field was taken to be the cosmic background at 3 K. We used a linear interpolator to estimate the line intensities for parameters in the intervals to better constrain the parameters and to allow for a continuous examination of parameter space. Nevertheless, the accuracy of the best-fit parameters remains limited by the resolution of the grid.

We fitted the A- and E-CH$_3$OH (5-4) ($v_t = 0$) line profiles with Gaussian models pixel-by-pixel for all the sources. We assumed that the linewidth is identical for all the K components

and for A and E types. We also took into consideration that the HNCO 11$^{0,11}$–10$^{0,10}$ line is blended with the CH$_3$OH 5$_{−}$1–4$_{−}$1 E line ($\delta V \sim 8.44$ km s$^{-1}$) and removed the contribution from this HNCO line by fitting an additional Gaussian component of which the linewidth and amplitude are free parameters.

We employ the Markov chain Monte Carlo (MCMC) method with an affine invariant sampling algorithm$^6$ (emcee, Foreman-Mackey et al. 2013) to perform the pixel-by-pixel fitting. Our prior assumption for each pixel adopted the temperature $T_{\text{rot}}$ predicted from the temperature profiles in Section 3.4, which means that for $T_{\text{kin}}$ we assumed a normal distribution centering at $T_{\text{rot}}$ and with a standard deviation of 5 K. For the other parameters, the priors were assumed to be uniform distributions. We used a likelihood distribution function which takes into account observational thresholds,

$$P \propto \prod_i p_i \prod_j p_j, \qquad (C.1)$$

where $p_i$ stands for probabilities of the ith data that is a robust detection and $p_j$ the jth data that gives a constraint by an upper limit. We adopt the normal distribution as likelihood function,

$$p_i \propto \exp \left[ -\frac{1}{2} \left( \frac{f_i^{\text{obs}} - f_i^{\text{model}}}{\sigma_i^{\text{obs}}} \right)^2 \right] \Delta f_i, \qquad (C.2)$$

$$p_j \propto \int_{-\infty}^{f_{\text{lim},j}^{\text{obs}}} \exp \left[ -\frac{1}{2} \left( \frac{f_j - f_j^{\text{model}}}{\sigma_j} \right)^2 \right] df_j, \qquad (C.3)$$

where $f_i^{\text{obs}}$ (or $f_j$) stands for the observed intensity (or intensity upper limit) obtained from Gaussian fit; $f_i^{\text{model}}$ (or $f_j^{\text{obs}}$) the model intensity; $\sigma_i^{\text{obs}}$ the standard error of the observed intensity, which was adopted as the fitted 1$\sigma$ error of the Gaussian fit, $\Delta f_i$ being the data offset from the true value of $f_i$; and $df_j$ the integrated flux probability to the detection threshold $f_{\text{lim},j}^{\text{obs}}$.

The starting points (initialization) for the chains were chosen to be the parameter set corresponding to a global $\chi^2$ minimum calculated between the grid models and the observed values. We also employed the burn-in phase in the MCMC chains and several resets of the starting-points to ensure the final chains are reasonably stable around the maximum of the probability density. An example posterior distribution of the fitted parameters is shown in Figure C.1. The obtained $n(H_2)$ maps and column density maps of CH$_3$OH for all sources are shown in Figure 9 and Figure C.2, respectively.

After the first run of generating parameter maps, we find that the $n(H_2)$ of clump G31 are truncated to ~10$^9$ cm$^{-3}$, the parameter boundary of our conducted RADEX models. Therefore we also ran a larger grid with $n(H_2)$ up to ~10$^{11}$ cm$^{-3}$ with the same range of $T_{\text{kin}}$ and $N_{\text{mol}}$ as the first grid, and re-derived the parameter maps for this clump.

## Appendix D: LTE analysis for other lines

### Appendix D.1: CS/C$^{34}$S (5-4)

CS transitions are commonly adopted tracers for dense gas. Our spectral setup covered the J=5-4 transition of CS and its rare isotopolog C$^{34}$S for the seven sources. For G10 the data we utilized only covered the C$^{34}$S (5-4) line. In the observed sources the CS

<footnote>$^6$ A detailed description of the method can be found in the emcee documentation.</footnote>
<table><caption>**Table B.1.** Information of the centimeter continuum data collected from NRAO archive.</caption><thead><tr><th>Source</th><th>Freq. (GHz)</th><th>Resolution (′′)</th><th>Obs. Date</th><th>Integrated flux (Jy)</th><th>Extrapolated total flux of free-free emission at 1.2 mm (Jy)</th></tr></thead><tbody><tr><td>G31</td><td>43.3</td><td>1.62</td><td>2001-11-26</td><td>0.42</td><td>0.40</td></tr><tr><td>G08b</td><td>8.5</td><td>0.85</td><td>2005-02-28</td><td>0.99</td><td>0.70</td></tr><tr><td>G10</td><td>23.8</td><td>0.09</td><td>2002-02-01</td><td>3.55</td><td>2.80</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Fig. C.2.** CH₃OH column density maps ($N_{CH_3OH}$) derived from RADEX modeling of all target sources. The beam of CH₃OH 5-1-4-1 E line is indicated in the bottom left corner. Gray contours indicate the SMA 1.2 mm continuum level from 0.1 to 0.9 × peak flux in five steps of uniform interval.</figcaption></figure>

line and the 1.2 mm continuum emission trace similar structures. Additionally, the CS line images revealed some spatially more extended structures.

Assuming that CS and C³⁴S trace gas of identical temperature and linewidths, the optical depth of the C³⁴S line can be estimated by

$$ \frac{W_{\text{CS}}}{W_{\text{C}^{34}\text{S}}} = \frac{[1 - \exp(-R\tau_{\text{C}^{34}\text{S}})]}{[1 - \exp(-\tau_{\text{C}^{34}\text{S}})]}, \quad (\text{D.1}) $$

where $W_{\text{CS}}$ and $W_{\text{C}^{34}\text{S}}$ are the integrated line intensities, and we assumed a ³²S/³⁴S abundance ratio of $R = 22.6$, which is consistent with terrestrial value (e.g., Lodders 2003), although there are indications that the ratio varies as a function of galactocentric distance (e.g., Humire et al. 2020). The CS line, as expected, is broader than the C³⁴S line, and shows self-absorption profiles in most of the cases. These are two effects that can mitigate each other in biasing the $\tau_{\text{C}^{34}\text{S}}$ estimate as our derivation was based on the integrated intensity ratios. We found that $\tau_{\text{C}^{34}\text{S}}$ ranges from
less than 0.1 to 0.25 in most of the sources, while it reaches 0.3-0.45 in G08b and values as high as 0.9 in G31.

Assuming a beam filling factor $f = 1$, a lower limit of the excitation temperature ($T_{ex}$) can be estimated from the brightness temperature of C$^{34}$S,

$$J(T_b) = [J(T_{ex}) - J(T_{bg})][1 - \exp(-\tau_{C^{34}S})]. \quad (D.2)$$

Since the $n_{crit}$ of CS is higher than the thermometers except CH$_3$CN used in Sect. 3.4, and because there is more extended cold envelope gas contributing to the emission of CS line, as expected, the derived $T_{ex}$ is in general lower than the $T_{rot}$ profiles derived from the observations of multiple transitions of the other molecular thermometers (see Sect. 3.4). The derived $T_{ex}$ is <10 K in G18, in the range of 20-30 K in most of the other sources, and reaches ~45 K at the emission peak of G31.

Under an LTE assumption, the overall column densities of the C$^{34}$S ($N_{C^{34}S}$) can then be estimated by evaluating the partition functions, adopting $T_{ex}$, with

$$N_{tot} = \frac{8\pi\nu^3}{c^3 A_{ul}} \frac{Q_{rot}}{g_u} \exp\left(\frac{E_{up}}{kT_{ex}}\right) \left[\exp\left(\frac{h\nu}{kT_{ex}}\right) - 1\right]^{-1} \int \tau_{C^{34}S} dv, \quad (D.3)$$

where $g_u$ stands for the level degeneracy and $Q_{rot}$ the partition function at $T_{ex}$. The derived C$^{34}$S column densities are in the range of $4.4 \times 10^{12} - 5.0 \times 10^{14} \text{ cm}^{-2}$.

In an alternative approach, we assumed that the $T_{ex}$ of C$^{34}$S is the same with the $T_{rot}$ measured from the other molecular thermometers, and then directly solved for $\tau_{C^{34}S}$ based on an LTE assumption. The $N_{C^{34}S}$ derived this way is 1.1–1.5 times larger for warmer sources except for G18. In G18 the $N_{C^{34}S}$ derived with this approach is ~10 times smaller than that derived based on the CS/C$^{34}$S intensity ratio. This is due to the dependence between $N_{C^{34}S}$ and $T_{ex}$ under LTE assumption: $N_{C^{34}S}$ drops significantly below ~30 K and slowly increases for larger values of $T_{ex}$.

## Appendix D.2: CCH (3-2)

The hyperfine line components of CCH resulting from electron-nucleus interactions allow for a direct measure of the CCH line optical depth. We followed a fitting procedure similar to that described in Estalella (2017), adopting relative line intensities from LTE predictions. From hyperfine line fitting we derived excitation temperature and optical depth on a pixel-by-pixel basis, and then used these values to estimate the CCH column densities. We found that the optical depth of the main line $\tau_m \sim 0.1$ in extended regions, but can become optically thick at high column densities in localized regions. In most of the sources, $\tau_m$ reaches ~3, but in G18 $\tau_m$ reaches ~10 at the center, which is close to the upper limit in our fitting procedure. The column density of CCH ranges between $1.2 \times 10^{14}$ and $6 \times 10^{15} \text{ cm}^{-2}$ across the emission region for all the sources. In general, the CCH column density distribution shows porosity over extended regions (~0.2-0.3 pc) except G18, and is systematically larger in the outer regions, resembling a ring-like structure, especially for clump G08b and G31. The central region of the hot massive core G31 is almost completely devoid of CCH emission. Jiang et al. (2015) presented high angular resolution CCH observations toward clump G10 together with three other more evolved high-mass clumps than those in our sample, which also exhibit ring-like distributions.

## Appendix D.3: SO and SO$_2$

Our spectral setup covered three SO lines (Table 3) which have similar upper level energies $E_{up}$. These three lines were

detected in all clumps except G18. The resolved SO intensities were spatially compact and were confined to the area where 1.2 mm continuum emission was detected. The SO$_2$ $14_{0,14}-13_{1,13}$ ($E_{up} \sim 94$ K) line was detected in the all clumps except G18 and G28; in G28 a lower excitation transition SO$_2$ $3_{2,2}-2_{1,1}$ ($E_{up} \sim 15$ K) was marginally detected.

These SO and SO$_2$ transitions are likely thermalized given that the emission closely follows the extension of the central core, and are of relatively high critical densities (Table 3). We cannot directly derive the excitation temperatures of the SO and SO$_2$ molecules because they do not cover multiple transitions, or because the covered transitions have similar $E_{up}$. Therefore, we adopted the gas temperature derived from the thermometer lines (Section 3.4) when estimating their column densities, assuming LTE conditions.

## Appendix E: Radiative transfer modeling of multi-wavelength continuum

We divided the 870 µm (using APEX/LABOCA), 450 µm (using JCMT/SCUBA-2), or 350 µm (using CSO/SHARCII or APEX/SABOCA) images (whenever available) of each clump into annuli that have intervals equal to half of the beam FWHM, and then (projected) radially averaged the intensities in the annuli. When the intensity distribution is largely asymmetric with respect to the clump center (e.g., when there is a bright adjacent core or clump or a bright external gas filament), we trimmed the contribution from those asymmetric (sub)structures before making the averages. In practice, at each radius we first derived the mean and standard deviation ($I_{std}$) of the intensities and then masked the pixels at which the sum of the radial intensities deviates from the mean by more than 0.9 times the sum of radial $I_{std}$.

We used radiative transfer models to invert these derived radial intensity profiles to radial density profiles. We employed the publicly available Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012). To find the simulated intensity profiles that match the observed ones, our modeling ergodically visited the parameter spaces of the assumed functional form for density (Equation 2, Sect. 3.6). This means that the resultant models, when their number is sufficiently large (we adopt 10,000 models), can represent the average statistical properties of the models constructed from the full parameter space.

To benchmark and refine the radial temperature profile, we conducted aperture photometry with multi-wavelength dust continuum data. We used the tools in Python package photutils (Bradley et al. 2021). The total fluxes at each wavelength of 24, 70, 160, 250, 350, and 500 µm from Spitzer/MIPS (Carey et al. 2009), Herschel/PACS, and Herschel/SPIRE images are measured and they provide an observed SED profile. We then let RADMC-3D generate the SEDs ranging from mid-IR to millimeter wavelengths and ensure that the best-fit model can reproduce an SED profile which is broadly consistent with the observed SED.

We assumed that dust opacity does not have spatial variation and quoted the opacity model from Ossenkopf & Henning (1994). Specifically, we quoted the column evaluated for the thin ice mantle coated dust which was coagulated for $10^5$ years in an environment with $10^6 \text{ cm}^{-3}$ gas density (hereafter OH5 model). This model has been successfully applied in previous studies to explain the radial profiles of dust emission around high-mass embedded protostars (Mueller et al. 2002; Rolffs et al. 2011). The maximum grain size in the OH5 opacity model is under the Rayleigh limit, and thus the scattering opacity can be neglected in the modeling.
The hydrogen gas density profile ($\rho_{\text{bulk}}(r)$) is described by the functional form of Equation 2. Using a Monte Carlo radiative transfer to evaluate the temperature distributions (e.g., based on assumptions of heating sources and interstellar radiation field) is subject to a very large degree of freedom. Instead, we adopted the radial gas temperature profile probed by the thermometers in Section 3.4, following the form of Equation 3.

The overall mass of individual clumps derived by single-component SED fitting ($M_{\text{sed}}$ as in Table 2) is used as a prior for our subsequent modeling. Given that we only need to fit two free parameters ($\bar{\rho}$, $q$, cf. Equation 2), the MCMC method has no clear advantage. Instead, for each observed molecular clump, we drew 10,000 samplers from the parameter space of $q = 0.0-2.5$, $\bar{\rho}=10^3-10^6$ cm$^{-3}$ with a random process. The likelihood function of the clump total mass follows the truncated normal distribution (e.g., Lynch 2007):

$$ \ln\mathcal{L} = \begin{cases} -0.5\left(\frac{M-M_{\text{sed}}}{3M_{\text{sed}}}\right)^2 & (0.2M_{\text{sed}} \le M \le 5M_{\text{sed}}) \\ -\ln f & (M < 0.2M_{\text{sed}} \text{ or } M > 5M_{\text{sed}}) \end{cases} \quad (\text{E.1}) $$

We then ran RADMC-3D for each sample and convolved the derived images with the corresponding Gaussian beams for comparison with the multi-wavelength observations. When summing the $\chi^2$, we assumed that the observational data can have nominal ~20% errors, and thus adopted one standard deviation in the radial profile calculation as the uncertainty of the observational data. The $\chi^2$ calculation follows

$$ \chi^2 = \sum_{\lambda} \sum_{r} \frac{(f_{\text{mod}} - I_{\text{obs}})^2}{I_{\text{std}}^2 + (0.2I_{\text{obs}})^2}, \quad (\text{E.2}) $$

where $\sum_{\lambda}$ and $\sum_{r}$ respectively denote summing over all wavelengths and sampled radii, and $f_{\text{mod}}$ stands for the model intensity at a certain radius.

Figure 10 compares the observed radial intensity profiles with the best-fit models. The fitted values of $q$ and $\bar{\rho}$, the predetermined clump radius, overall clump mass, and the characteristic density determined from the best-fit model at 0.1 pc radius are summarized in Table 6. The posterior probability distributions of $\bar{\rho}-q$ are shown in Figure E.1, with the best-fit parameter set indicated. In Figure 11 the observed multi-wavelength flux densities are compared with the SED generated from best-fit RADMC-3D models.

It can be seen that the observed SEDs of sources G18, G19, G08a, G13, and G31 show large deviations from those of the best-fit models that were produced based on the assumed $T(r)$ profiles. For sources G18, G19, G08a, and G13 the deviations are partially expected as these three sources exhibit the lowest densities in the sample, from the clump center to the outer region. It is likely that $T_{\text{rot}}$ estimated from the aforementioned thermometers is biased to the small proportion of the dense gas at each radius (clumpiness), while the temperature for the bulk gas at each layer or the average gas temperature is lower. In addition, for the intermediate-scale gas having densities $\le 10^5$ cm$^{-3}$, which are mostly probed by CH$_3$CCH (12-11) and H$_2$CS lines in terms of gas temperature, the higher gas temperature than dust temperature may also have an origin from turbulent heating, as under such gas densities thermal coupling between dust and gas is weaker (Pan & Padoan 2009). For G31, which is much denser in terms of bulk gas radially but shows the most prominent monolithic core in the center, the overestimation of temperature likely originates from the optical depth effect.

To refine $T(r)$ as guided by the observed SED shape, we retain the parametric form (Equation 3) as elaborated in Section 3.6, and only adjust the parameter $r_{\text{in}}$ (Equation 3), scaled by a factor <1. Based on the adjustment, we regenerate SED profiles from RADMC modeling and find the best-fit rescaled $r_{\text{in}}$, which is listed in Table 6 for the five sources. For these five sources we then iterate the fitting of radial density profiles from RADMC-3D calculations, based on the adjusted temperature profiles $T(r)$. From Figure E.1, the best-fit parameter set before and after adjusting $T(r)$ are shown together. Decreasing $r_{\text{in}}$ in the temperature form is equivalent to reducing the steepness and absolute value of temperature radial profile, which results in an increment of mean gas density and density profile slope in the radial intensity profile fits, as expected. As a further benchmark, we use RADMC-3D to self-consistently calculate the dust temperature. To convert the clump luminosity to a central stellar source, we use the stellar evolution models of solar metallicity from Choi et al. (2016), which give relations between luminosity, mass, radius, and temperature for ZAMS to estimate the stellar $T_{\text{eff}}$. The re-iterated best-fit density model is used to describe the envelope structure. The resultant SED is also shown in Figure 11.

## Appendix F: Radiative transfer modeling of CH$_3$OH and CH$_3$CCH lines: Benchmarking the results from one-component non-LTE and LTE models

The RADEX analysis of CH$_3$OH (5-4) lines in Section 3.5, with one-component non-LTE assumption, can be biased due to the mixed contribution along the LOS. On the other hand, the RADMC-3D modelings for dust continuum in Appendix E better represents the enclosed masses within certain radii by constraining $\rho_{\text{bulk}}$, while they cannot provide constraints on the higher volume density internal structures. To refine our estimates of the physical parameters, we built on these two efforts to conduct 3D radiation transfer forward modelings of CH$_3$OH (5-4) and CH$_3$CCH (12-11) lines. Specifically, we performed the non-LTE modeling for the excitation conditions using the LIME code (Brinch & Hogerheijde 2010) and then compared the synthetic spectra with the observed spectral cubes. The spatial distributions of CH$_3$OH (5-4) and CH$_3$CCH (12-11) line emission are rather extended such that they better characterize the majority of dense gas in the clumps. We focus on benchmarking the radial density profile of the dense gas ($\rho_{\text{dense}}$). We fixed the gas temperature profiles ($T(r)$; Section 3.6) according to the profiles obtained in Section 3.4 (and refined in Appendix E), when constructing the input models. As $T(r)$ has been verified by full radiative transfer of dust continuum by building SEDs, the CH$_3$CCH (12-11) modeling conducted here is further used as a sanity check that the assumed temperature profile can reproduce the line emission of this thermometer.

The gas density radial profiles were initially specified to be the best fits from the modeling continuum intensity profiles, $\rho_{\text{bulk}}$ (Equation 2, see details in Appendix E). Based on this density model and temperature profile $T(r)$, we then experimented with various assumptions of the abundance spatial variations of CH$_3$OH and CH$_3$CCH, motivated by the results of the chemical network from Belloche et al. (2017) and Calcutt et al. (2019).

Based on the chemical network of Belloche et al. (2017), Calcutt et al. (2019) calculated the variation in gas-phase abundance of CH$_3$CCH as a function of warm-up time at different densities which range from $10^7$ to $10^{10}$ cm$^{-3}$. In these calculations, a two-stage physical evolution is assumed: a cold collapse stage is followed by a static warm-up stage that reaches a