Cloud Structure of Three Galactic Infrared Dark Star-forming Regions from Combining Ground- and Space-based Bolometric Observations

High angular resolution, ground-based continuum observations at 350 μm toward the three IRDCs G11.11−0.12, G28.34+0.06, and G14.225−0.506 were carried out using the SHARC2 bolometer array, installed on the CSO Telescope (PI: H. B. Liu). The array consists of 12 × 32 pixels.¹³ The simultaneous field of view (FOV) provided by this array is 259 × 097, and the diffraction limited beam size is ∼88. The data of G11.11−0.12 were acquired on 2014 March 28 (${\tau }_{225\mathrm{GHz}}$ ∼ 0.06), with an on-source exposure time of 90 minutes. G14.225−0.506 and G28.34+0.06 were observed on 2014 March 27 (${\tau }_{225\mathrm{GHz}}$ ∼ 0.05), with 80 and 50 minutes of on-source exposure time, respectively. The telescope pointing and focusing were checked every 1.5–2.5 hr. Mars was observed for absolute flux calibration. We used the standard 10' × 10' on-the-fly box scanning pattern, and the scanning center for each source is listed in Table 1. Basic data calibration was carried out using the CRUSH software package (Kovács 2008). We used the -faint option of the CRUSH software package during data reduction, which optimized the reconstruction of the faint and compact sources with the cost of the more aggressive filtering of extended emission. Nevertheless, the extended emission components will ultimately be complemented by the observations of space telescopes (see Section 2.4 for more details). The final calibrated map was smoothed with a Gaussian kernel with 2/3 beam FWHM (-faint option) to an angular resolution of 96 for an optimized sensitivity and source reconstruction, and the rms noise levels we measured from the approximately emission-free areas of 56, 53, and 60 mJy beam⁻¹ for G11.11−0.12, G28.34+0.06, and G14.225−0.506, respectively.

The 870 μm submillimeter continuum observations toward G14.225−0.506 were conducted (PI: G. Busquet) using the LABOCA bolometer array, installed on the Atacama Pathfinder EXperiment (APEX).¹⁴ The array consists of 259 channels arranged in nine concentric hexagons around the central channel. The FOV of the array is $11\buildrel{\,\prime}\over{.} 4$, and the angular resolution of each beam is $18\buildrel{\prime\prime}\over{.} 6\pm 1^{\prime\prime}$. The observations were carried out on 2008 August 24 and 31 under weather conditions with zenith opacity values that ranged from 0.15 to 0.24 at 870 μm. The observations were performed using a spiral raster mode mapping, providing a fully sampled and homogeneously covered map in an area of $15^{\prime} \times 15^{\prime}$.

Calibration was done using observations of Mars as well as secondary calibrators. The absolute flux uncertainty is estimated to be ∼8%. Pointing was checked every hour, finding an rms pointing accuracy of $2^{\prime\prime}$, and focus settings were performed once per night and during the sunset. The data was reduced using the MiniCRUSH software package (see Kóvacs 2008). The data reduction process included flat-fielding, opacity correction, calibration, correlated noise removal, and de-spiking. Further details of the data acquisition and data reduction are described in Busquet et al. (2016).

We retrieved the available night observations¹⁵ of JCMT Submillimetre Common-user Bolometer Array 2 (SCUBA2; Chapin et al. 2013; Dempsey et al. 2013; Holland et al. 2013) at 850 μm from the online data archive. (Program ID is M11BEC30 for G11.11−0.12 and G28.34+0.06.) Ancillary data also includes level 2.5 and level 3 processed, archival Herschel¹⁶ images, which were taken by the Herschel Infrared Galactic Plane (Hi-GAL) survey (Molinari et al. 2010) at 70/160 μm using the PACS instrument (Poglitsch et al. 2010) and at 250/350/500 μm using the SPIRE instrument (Griffin et al. 2010). Observation IDs for G14.225−0.506 are 1342218997 and 1342219000, for G11.11−0.12 they are 1342204952 and 1342218965, and for G28.34+0.06 it is 1342218694. Planck/High Frequency Instrument (HFI) 353 GHz images are also used. For our combination purpose, we use the Herschel SPIRE extended emission products from the data base. The Planck 353 GHz images, which are in units of ${{\rm{K}}}_{\mathrm{CMB}}$, are converted to Jy beam⁻¹ (Planck HFI Core Team et al. 2011a, 2011b; Zacchei et al. 2011).

Our procedure to combine images taken from ground-based and space telescope observations in the Fourier domain, and then iteratively derive the high angular resolution dust temperature and dust/gas column density images, is similar to what was introduced in Lin et al. (2016). Since the observed IRDCs for the present paper have relatively low temperatures, it is particularly important to precisely determine the long wavelength part of the spectrum. Therefore, we additionally included the following steps to improve the quality of the combined 850 μm image before making our final modified blackbody SED fits.

We extrapolated the 850 μm flux from SED fits to Herschel PACS 160 μm and SPIRE 250/350/500 μm. In these fits, we fixed the dust emissivity index β to 1.8, which is the mean value measured for the Galactic disk (Planck Collaboration et al. 2011, 2014). The Herschel-extrapolated 850 μm image has the same resolution as SPIRE 500 μm, ∼37''. Then we proceeded to combine the above mentioned 850 μm image with Planck 353 GHz, and then further combined the result with the SCUBA2 850 μm image. In this way, the >5' scale and the <3'–4' scale structures are dominantly constrained by the Planck and the SCUBA2 images, respectively. The role of the Herschel-extrapolated 850 μm image is to complement the small range of spatial scales, which is poorly sampled by Planck and SCUBA2.

Figures 1–3 show the derived dust temperature and gas column density maps, assuming a gas to dust mass ratio of 100. For detailed calculations and procedure, see Appendix A and Lin et al. (2016). The derived gas masses of these three IRDCs range in ∼3–9 × 10⁴ ${M}_{\odot }$. The overall gas masses and bolometric luminosities derived from these maps are summarized in Table 1 .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The column density maps of these IRDCs present considerably different overall morphologies. On the ∼30 pc scale, IRDC G11.11−0.12 appears as an integral shaped filament, which shows some wiggles on ∼10 pc scales. In addition, the high angular resolution we achieve reveals that the large-scale filament may in fact consist of bundles of filaments. The majority of dense gas structures in G11.11−0.12 have dust temperature lower than 20 K. We are able to identify one internally heated source embedded approximately at the center of the overall filament, where previous observations have found water masers, Class II methanol masers, and compact centimeter continuum sources (Pillai et al. 2006; Rosero et al. 2014; Wang et al. 2014). There are several isolated heating sources, which are projected adjacent or well outside of the dense filament. How many of those heating sources are directly associated with this molecular cloud is not yet certain.

IRDC G14.225−0.506 is resolved into a web of dense gas filaments, which seem to align in two preferred directions (see also the discussion in Busquet et al. 2013). The two most significant column density peaks in Figure 2 correspond to the two previously known relatively massive (∼10³ ${M}_{\odot }$), ∼0.5 pc scale molecular clumps.¹⁷ The northern massive molecular clump, namely Hub-N, is adjacent to a compact H ii region, which presents a high dust temperature over 1'–2' angular scales. The rest of the dense structures have ∼18–20 K dust temperature in general. There are several heated sources embedded in the filamentary structures in G14.225−0.506.

The dominant dense gas structures in G28.34+0.06 can be described by a massive molecular clump in the north (P2, see Carey et al. 2000; Wang et al. 2008; Zhang et al. 2009), which is connected with a dense gas filament in the south. The southern gas filament has several embedded clumps.

Star-formation in the P2 clump is already active. The P2 clump is internally heated by the embedded young stars, and shows higher than 25 K dust temperature. Dust temperature in the southern molecular clumps remains under 20 K.

To systematically quantify the dense gas distributions in the observed IRDCs, and thereby permit comparisons with other observations and theoretical models, we perform analyses of the column density probability distribution functions (N-PDF), column density complementary cumulative distribution function (N-CCDF), and the two-point correlation functions (2PT) of gas column density.

The N-PDFs of the three IRDCs are presented in Figure 4. For the sake of enabling quantitative comparisons with other observations, we approximate the components in the N-PDFs that have a decreasing slope at lower column densities by a log-normal distribution function, and approximate the components that have a near constant slope at higher column densities by a power-law distribution function. In the later discussion, we will refer to these as the "log-normal" and "power-law" components. We emphasize that the functional forms we selected are merely approximations for the actual observed N-PDFs, which instead show richer features that cannot be described by a simple fit to these functional forms. Linking these functional forms to the underlying physics is not trivial, and one should not over-emphasize on the exact functional forms. We elaborate further on this point in Section 4.2 by examining comparisons of numerical hydrodynamic simulations.

Zoom In Zoom Out Reset image size

We fit an overall power-law and a power-law starting from the optimal column density cutoff based on Maximum Likelihood Estimation (MLE; Clauset et al. 2009; Alstott et al. 2014) to the N-PDF of G28.34+0.06. The results of these two fits appear consistent. At the high column density end, however, the N-PDF of G28.34+0.06 shows an excess over the fitted power law. The N-PDFs of G14.225−0.506 and G11.11−0.12 are better described by a log-normal component in the lower column density part, in addition to the optimal power-law fit at the high column density end. Whether the log-normal component is included in the fit or not does not affect the results of the optimal power-law fit because of the MLE approach used to derive the optimal cutoff for power-law fit instead of jointly fitting the two functional forms.

The N-PDF of G14.225−0.506 presents an apparent deviation from the log-normal and power-law fits over a wide range of the normalized column density at η ∼ 1–3, and therefore it is hard to define whether there is an excess or deficit at any range of η. In particular, a significant "bump" is seen around η ∼ 2.0–2.3. The N-PDF of G11.11−0.12 decreases the most rapidly with increasing η. The fitted parameters are tabulated in Table 2 and are also labeled in each N-PDF panel in Figure 4.

Table 2.  Fitting Results of Column Density Probability Distribution Functions (N-PDFs)

Target Source	Measured Threshold	Mean Column Density	Log-normal and Power-law Fits
	log10(N(H2))	log10(N(H2))	${\sigma }_{\eta }$	μ	s0	s1
G28.34 + 0.06	21.70	21.73	⋯	⋯	−3.98(0.01)	−3.86(0.02)
G14.225−0.506	21.80	22.04	0.43	21.99	⋯	−4.12(0.01)
G11.11−0.12	21.95	22.09	0.26	22.07	⋯	−6.75(0.09)

Download table as:  ASCIITypeset image

The two-point correlation (2PT) function we used follows the same form as that in Kleiner & Dickman (1984), but instead of calculating the correlation of column density fluctuations, we directly measure the correlation of column densities for each source across the observed field. The correlation strength at a separation scale (lag) of l is calculated by

$$\begin{eqnarray}&&{S}_{\mathrm{tr}}(l)=\displaystyle \frac{{\langle X({\boldsymbol{r}})X({\boldsymbol{r}}+{\rm{l}})\rangle }_{{\boldsymbol{l}}}}{{\langle X({\boldsymbol{r}})X({\boldsymbol{r}})\rangle }_{{\boldsymbol{l}}}},\end{eqnarray} \tag{ 1 }$$

where $X({\boldsymbol{r}})$ denotes the column density value at position r, and the angle brackets are the average over all pairs of positions with a separation of l. The final form of the correlation function is normalized by the peak correlation strength to enable comparisons between different fields. For a more detailed description, we refer to Lin et al. (2016).

Figure 5 shows the 2PT functions of the three IRDCs overplotted with those of the active OB-cluster-forming molecular clouds taken from Lin et al. (2016). The 2PT functions of the three IRDCs, in general, show a smooth decay of correlation strengths over all spatial scales. This is in contrast with most of the active OB-cluster-forming regions, which rapidly decrease at small scales. The dominant elongated cylindrical filaments in G11.11−0.12 and G14.225−0.506 may naturally have flattened 2PT functions since, by definition, filaments have comparable column density distributions over a wide range of separations. For G28.34+0.06, we spatially resolved complex local over-densities scattered across the column density map. The multiple fluffy and scattered over-dense substructures of this molecular cloud more closely resemble the cloud morphology of W43-south and G10.2−0.3, which also have similar 2PT functions (see Lin et al. 2016).

Zoom In Zoom Out Reset image size

The homogeneity of the column density fields could be estimated by the range of plateaus in the 2PT functions. In light of this, G11.11−0.12 would present the least mass concentration since its 2PT function remains at an almost constant correlation strength of up to ∼6 pc. G14.225−0.506, which harbours two pronounced massive clumps, shows a steep decrease at lags of less than 1 pc, indicating a relatively significant concentrated distribution at small scales. The 2PT function of G28.34+0.06 falls between those of G11.11−0.12 and G14.225−0.506 with a slight decay at small scales ($\lt 1\,\mathrm{pc}$) and a constant correlation plateau to ∼4 pc. We note that the 2PT function plateaus on the few parsec scales may be related to the elongated nature of cloud structures.

Figure 6 presents the column density complementary cumulative distribution function (N-CCDF) of these molecular clouds (this is also defined as the dense gas mass function in some previous works; see details in Kainulainen & Tan 2013; Ginsburg et al. 2015). We exclude W49A, in spite of its highest fraction of dense gas, to avoid a bias from its larger distance (∼11 kpc). The dense gas fraction is closely associated with the star-formation rate, raised based on observations toward both galactic and extragalactic star-forming regions (Gao & Solomon 2004; Wu et al. 2005; Lada et al. 2010). In Figure 6, we plot the N-CCDFs for all of the sources. The dashed gray reference lines indicate exponential decay of indexes between −0.06 to −0.12. The N-CCDF of OB-cluster-forming molecular clouds within our sample of IRDCs are generally more flattened in the high column density ends, except for G10.2−0.3 and G10.3-0.1. These sources hold a considerable mass fraction in the high column density regime. Despite its relatively high total mass, G11.11−0.12 has a sharp decrease of very dense gas and is consistent with its steep N-PDF power-law tail. We additionally plot two reference vertical lines representing a 10⁵ cm⁻², ∼0.1 pc core and a 10⁶ cm⁻², ∼0.03 pc core lying at the same distance.

Zoom In Zoom Out Reset image size

The N-CCDF result we obtained for G11.11−0.12 is consistent with Herschel results obtained by Kainulainen et al. (2013b), where the authors also use near- and mid-infrared absorption to derive an $\sim 2^{\prime\prime}$ column density map of this cloud. The large discrepancy of our column density map with their extinction derived one may mainly come from the resolution difference (and also the reprojection to a larger distance in our case). An ∼0.035 pc resolution extinction map enables detections of embedded dense cores, which are expected to give a significant rise to N-CCDF in the high column density end. We also notice the fact that the near- and mid-infrared extinction method is limited to detecting a lower range of column density (${A}_{v}\sim 100$, according to Kainulainen & Tan 2013). Therefore, to evaluate the origin of this difference in the high column density regime, one should perform systematic error estimates of the extinction method, which is beyond the scope of our current analysis.

The sharp decrease for G28.34+0.06 and G14.225−0.506 at less than ${\rm{N}}({{\rm{H}}}_{2})\sim 4.0\times {10}^{22}\,{\mathrm{cm}}^{-2}$ and the near flat region up to ${\rm{N}}({{\rm{H}}}_{2})\lt \sim 8.0\times {10}^{22}\,{\mathrm{cm}}^{-2}$ may place these clouds at an evolutionary stage in between G11.11−0.12 and OB-cluster-forming samples. These results imply that these two IRDCs already have a significant concentration in mass at relatively high column densities.

In this section, we compare the three observed IRDCs in terms of their N-PDFs results, and how these are linked to possible physical processes at work. We also compare their measured properties with observations of seven more evolved OB-cluster-forming regions reported in Lin et al. (2016). The masses and luminosities of all the sources are summarized in Figure 7. We caution that small sample size limits our analysis to being conjecturing rather than conclusive.

Zoom In Zoom Out Reset image size

The N-PDFs of G11.11−0.12 and G14.225−0.506 have a more prominent log-normal component than OB-cluster-forming regions. Despite the possible bias induced by non-uniform sensitivity and different distances toward different sources, the N-PDF results in this paper and in Lin et al. (2016) are derived only for values greater than the measured column density threshold of each region and normalized to enable reliability in the comparisons. The log-normal plus power-law tail N-PDFs of G14.225−0.506 and G11.11−0.12 and their highly filamentary structures resemble, but are not uniquely explained by, the self-gravitating cylinder model (Myers 2015). The relative significance of the high-density part in the N-PDF of G14.225−0.506 is higher than in G11.11−0.12, with a larger ${\sigma }_{\eta }$ and a shallower power-law tail index α. The increased spatial density concentration is also seen in the sources' 2PT functions: G14.225−0.506 has a more prominent decrease of correlation strength at small spatial scales. Turbulent flows or shocks in the large-scale, low column density regime may shape the initial morphology of these clouds, producing the log-normal component in the N-PDFs (Vazquez-Semadeni 1994; Klessen 2000). For example, Tackenberg et al. (2014) suggested that a large-scale accretion flow exists along the integral filament of G11.11−0.12. Unlike the extreme cases of W49A and G10.6−0.4, which may be undergoing global and local collapse simultaneously in the clouds, the development of power-law tails at the high column density ends of G11.11−0.12 and G14.225−0.506 may suggest that they are only gravitationally unstable over a small fraction of the observed area.

We note that the N-PDFs of these two IRDCs may not necessarily be simply understood as being caused fully by the interplay between turbulence and self-gravity, considering that they both have been suggested to have dynamically non-negligible magnetic fields (Pillai et al. 2015; Santos et al. 2016). The index of the N-PDF power-law tail and variance could be influenced by the magnetic field, as suggested by gravoturbulence simulations (Burkhart et al. 2015). In some of these simulations, the narrower standard deviation σ of the log-normal part of the N-PDF is found to be caused by magnetic fields acting as a cushion that sets column densities closer to their average value (Nakamura & Li 2008; Molina et al. 2012). The smaller variance and steep power-law tail of G11.11−0.12 seems to be compatible with these simulations. It has been suggested that the low column density end of the log-normal component may be biased due to limited image size (e.g., Lombardi et al. 2015). Our measurements might be affected by the completeness issue, hence, we avoid further discussions of the log-normal component in the present research.

The relation between changes in the N-PDF power-law tails and the evolutionary stage of molecular clouds is suggested by several recent observational and numerical simulation works. For example, Stutz & Kainulainen (2015) investigate the variations in the N-PDF slopes with young star-formation content in the Orion molecular cloud, and they find that the N-PDF slope is steeper in regions where there is a smaller fraction of Class 0 protostars. As the power-law tail flattens, the fraction of young protostars increases. The shallower power-law tail related to more active star formation is also suggested by Lombardi et al. (2015) based on investigations toward several molecular clouds. The possible relation between the power-law tail slopes and collapsing state of molecular clouds is also raised by Kritsuk et al. (2011), Ballesteros-Paredes et al. (2011), Burkhart et al. (2015) etc. in their simulations. Regarding cloud properties reflected by N-CCDF measures, Kainulainen et al. (2013a) find that the slope of N-CCDF relies sensitively on the turbulence driving mechanism and SFE, and magnetized simulations of clouds have steeper slopes than non-magnetized ones. The 2PT functions of the selected IRDCs show slowly decaying correlation strengths at small spatial scales, indicating less matter concentration in local areas compared with most of the more luminous OB-cluster-forming regions. On the other hand, the IRDCs exhibit a larger decrease of correlation strengths at larger scales. The different (mixed) driving modes and magnetized nature of the turbulent motions in these IRDCs, compared to those in luminous OB-cluster-forming regions, may be the origin of the above mentioned differences in our statistical measurements (Federrath et al. 2009).

The N-PDFs of G11.11−0.12 and G28.34+0.06 were also derived by Schneider et al. (2015) based on Herschel data only. Despite their different angular resolution, the overall shape of N-PDFs are similar, while we are able to resolve more localized dense structures. We note that the main caveat of utilizing N-PDF features to quantify cloud properties is that it is an indirect indicator of cloud physical state, unlike, i.e., volume density profiles. It is possible to convert column density maps to densities using simple assumptions of cloud geometry (Stutz & Gould 2016), but this is not readily applicable to our sources considering their complex morphological structures based on previous spectroscopic works, which show that they have multiple velocity components.

In Figure 7, we plot the total luminosity versus mass of these sources with reference lines corresponding to constant luminosity-to-mass ratios. The luminosity-to-mass ratio provides a measure of evolutionary state. The ratio places G28.34+0.06 in a similar stage with some of the OB-cluster-forming regions, whereas G14.225−0.506 and G11.11−0.12 are more quiescent regions. The slope of the N-PDF and the luminosity-to-mass ratio of the compared sources are presented in Figure 8. The data points in the plot suggest a correlation between these two quantities.

Zoom In Zoom Out Reset image size

We here compare our observations with recent numerical simulations by Dale et al. (2015), who followed the formation of OB clusters from molecular clouds, and examined the impact of the clusters' feedback on the clouds. The simulations were performed in the Smoothed Particle Hydrodynamics formalism. Smooth spherical clouds of a range of masses were initialized with a turbulent velocity field and allowed to evolve and form stars. Once a few massive stars had formed in each cloud, the effects of their winds and ionizing radiation were modeled. The evolution of the clouds followed as close as was practicable to the epoch when the massive stars were expected to explode as supernovae, allowing the influence of the winds and radiation to be isolated.

We choose the 10⁴ M_⊙ Run I calculation from Dale et al, since it is among those with the best linear resolution. We compute simple column-density images along the simulation z-axis of 10 × 10 pc subfields of the cloud, comparable to the physical sizes of our observed IRDCs (see also Appendix B). Column-density images are produced at four epochs prior to the formation of stars and the initiation of feedback (top row of Figure 9), at the point when feedback begins to act (first panel in the bottom row of Figure 9), and for three later epochs (remaining panels in the bottom row of Figure 9).

Zoom In Zoom Out Reset image size

The four images of the top row in Figure 9 show the development of self-gravitating filamentary structures from the cloud's turbulent velocity field before the onset of star formation. The gaseous filaments are initially not gravitationally unstable and, at early stages, their main role is to feed gas to a forming, parsec-scale massive molecular clump (the bright compact object located in the lower left corner of the two rightmost images in the top row). This quickly accumulates a large amount of mass and becomes the earliest site of star formation. In Figure 10, we compute N-PDFs for these column-density images in the same way as we did for the observational data. The absence of many massive molecular clumps in the simulated cloud is due to the initial dominance of global collapse.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At the earliest epoch when the filamentary structure is least well-defined, the PDF has a form close to log-normal in appearance. We caution that the drop-off at low η is likely an incompleteness effect caused by the limited FOV and the lowest density contours not being closed. As the simulation progresses and the filaments become self-gravitating in preparation for forming stars, the N-PDFs develop pronounced power-law tails, which extend to higher and higher values of η, as do the observed plots (Figure 10). The slopes of the power laws are initially very steep, but flatten with time, approaching the slopes in the observed N-PDFs.

In even more striking resemblance to the observations of G14.225−0.506 (Figure 4), the simulated N-PDFs also shows a similar bump in the simulated data at $\eta \approx 2\mbox{--}2.5$ ((b)–(d) in Figure 9). This is dominated by the protocluster core mentioned earlier, indicating that material piles up in this structure for some time, slowing its progression to still higher densities and resulting in a bump in the N-PDF. There is also a minor contribution to the bump from the ambient dense gas filaments, which are immediately feeding the massive molecular clump. The massive molecular clump rapidly contracts due to self-gravity, which produces internal structures with high density and high column density. It then no longer appears as a bump in the N-PDF, whose slope extends further to higher values ((e) and (f) in Figure 9). The N-PDF evolves to a broken power-law shape at larger than zero η, showing an excess of high column density pixels. Such an N-PDF is similar to that observed in G28.34+0.06, and those of some active OB-cluster-forming regions (e.g., G10.6−0.4, W49A) reported in Lin et al. (2016). We note that an intermediate density bump present at volume density probability function (ρ-PDF) in the cloud's early evolved phase is also seen in the formation of protoclusters simulated by Lee & Hennebelle (2016).

In our simulations, the log-normal like N-PDF is thus not smoothly lifted to become a power-law like N-PDF, but instead changes via a rather dynamical process that produces non-smooth and non-steady features. The bottom panels (e)–(h) of Figure 9 show the reaction of the N-PDF to feedback from the OB cluster, which eventually forms in the dense core. The massive molecular clump forms a condensed cluster of stars; some stars are ejected due to many-body gravitational interactions. Subsequently, the massive clump is destroyed by feedback and by consumption due to stellar accretion. Feedback of the stars formed in the massive clump then quickly disperse the low density gas and destroys the filamentary structure from which the clump originally formed. As a result, a cavity is gradually created at the center of the system. The action of feedback is initially rather modest and is essentially to drive the high-density tail to higher and higher values of η. The N-PDF at this stage shows a deficit of high column density gas when compared with a single power-law fit ((g) in Figure 9), which is qualitatively similar to the observed N-PDF of the evolved OB-cluster-forming region G10.2−0.3 (Lin et al. 2016). However, G10.2−0.3 presents several H ii bubbles, which may shape the N-PDF differently. Eventually, the power-law slope of high column density gas becomes steeper again ((h) in Figure 9), as feedback becomes intense enough to significantly shape the high column density end. The smoothed N-PDFs shown in Figure 10 are simulated images smoothed to the resolution of our achieved resolution $10^{\prime\prime}$ at a median distance of our sources of 3 kpc. The smoothed N-PDFs and their fitted optimal power-law tail slopes are close to the original ones except that they exhibit clear truncations, which means that smoothing and differing distances do not affect our conclusions significantly but that high-resolution column density maps are essential to robustly recover the suppressed power-law tails at the high column density end.

As discussed in detail in Dale et al, feedback in this cloud is able to generate dense gas by compressing the outer regions of the cloud into shells, but much of this gas fails to become gravitationally unstable because it is located far from the cloud's main potential well, and is stabilized by geometrical stretching and turbulence. We also examined one of the larger clouds modeled by Dale et al. (2015), the 3 × 10⁵ ${M}_{\odot }$ UUCL cloud). We obtained qualitatively similar results to those obtained for Run I, though the coarser mass and length resolution of Run UUCL does not allow the high-density end of the PDFs to be resolved as is possible in Run I.

In the field of high-mass star formation, the evolutionary stages of massive (e.g., ${M}_{\mathrm{gas}}$ > 10⁴ ${M}_{\odot }$) molecular clouds are commonly separated into the infrared dark stage, the active star-forming stage, and the dispersed stage. IRDCs are considered to be the youngest evolutionary stage, representing the initial conditions of high-mass star formation. However, unlike low-mass stars, high-mass stars form in much shorter (Kelvin–Helmholtz cooling) timescales than the characteristic (free-fall) timescales for molecular cloud evolution. Therefore, characterizing molecular cloud evolution based on the illumination from embedded high-mass stars, may be coarse and biased. In terms of terminology, whether a molecular cloud is infrared dark or not may also depend on the sensitivity and the FOV of the observations, and the time variability of the observed stellar emission, which can be ambiguous (e.g., Feng et al. 2016).

Establishing evolutionary tracks of massive molecular clouds based on their structures may be more essential and crucial for facilitating studies of the co-evolution of the (proto)stellar and core mass functions with the molecular cloud (e.g., Zhang et al. 2015). Based on a limited sample (the present work and Lin et al. 2016), and also based on the results of hydrodynamic simulations (Section 4.2), we hypothesize four common characteristic evolutionary stages of massive molecular clouds. From the earlier to the later evolutionary stages, they are (1) the cloud integration stage, (2) the stellar assembly stage, (3) the cloud pre-dispersal stage, and (4) the dispersed cloud stage.

In the cloud integration stage, high column density (e.g., ${\rm{N}}({{\rm{H}}}_{2})$ ≥ 10²² cm⁻²) gas structures are forming on >10 pc scales. At this stage, the N-CCDF remains steep. The gas/dust temperature may be <20 K and the gas temperature is, in general, anti-correlated with the gas/dust column density, despite a few hot spots probably associated with deeply embedded intermediate-/high-mass (proto)stars. Cloud structures at this stage are characterized by a steep and decreasing slope of the power-law tail in the N-PDFs, and the diffuse nature of the 2PTs do not show a significantly distinguishable component in the short lags. In ∼10 pc scales, the low column density part of the N-PDFs may be described by log-normal distributions, until the self-gravitational contraction outweigh other physical mechanisms present from the initial conditions (e.g., turbulence, magnetic field, etc.). In addition, the N-CCDF and the luminosity-to-mass ratios are both low.

In the stellar assembly stage, the slopes of the power-law tail in the N-PDFs become shallower due to the overall self-gravitational contraction of the molecular cloud. In addition, the highest column density end of the N-PDF may start to present an excess of single power-law tail due to the presence of very massive molecular clumps/cores, which are forming the highest mass stars of the final clusters. The 2PT functions at this stage start to present a distinguishable component at lags $\lesssim 1\,\mathrm{pc}$, which is related to the presence of the above mentioned very massive molecular clumps/cores. The luminosity-to-mass ratio at this stage increases very rapidly due to the quickly forming high-mass stars. From the two-dimensional histogram of the gas/dust temperature and column density, we can also identify positive correlations from the core/clump material of the immediate surroundings of the newly formed high-mass stars.

In the cloud pre-dispersal stage, the remaining dense gas associated with massive molecular clumps/cores is either accreted onto the young massive star(s) (Ginsburg et al. 2016), or dispersed by the stellar feedback (Dale et al. 2015). The distinguishable component at the $\lesssim 1\,\mathrm{pc}$ lags of the 2PT functions become less prominent or it fully disappears. The slope of the power-law tail in the N-PDFs remains shallower than that in the cloud integration stage. The luminosity-to-mass ratio and the dense gas fraction at this stage will remain high as well. It is until the dispersed cloud stage that radiative feedback and the expansion of the ionized gas become capable of reverting the contracting motions on all spatial scales, significantly re-shaping molecular gas structures. The N-PDF becomes steeper back again and may show an increasing slope toward the high column density end. This evolutionary stage can be easily distinguished from the cloud integration stage and other evolutionary stages with similar N-PDF slopes because the gas/dust temperatures are much higher, and the presence of extended mid-infrared and radio continuum emission, which trace the extended and expanding H ii regions. Whether the N-PDFs can be reverted by the feedback or not, may also depend on the initial cloud morphology (e.g., the cross-section against feedback). In the case in which the stellar feedback is only strong enough to disperse low density gas, the piled up high-density gas may re-collapse in some free-fall timescales, recycling again.

We tentatively suggest that the cloud integration stage and the dispersed cloud stage can be separated from the stellar assembly stage and cloud pre-dispersal stage by the s₁ = −4 power-law index of the N-PDF, and by the correlation strength at the short-lag break in the 2PT function of correlation strength 0.8. The cloud integration stage may be further separated from the three more evolved stages by a luminosity-to-mass ratio ∼40 ${L}_{\odot }$/${M}_{\odot }$. The stellar assembly stage may have a more prominent decrease of 2PT at small scales than the cloud pre-dispersal stage, with the latter possibly exhibiting a deficit of the singular power-law tail of N-PDF at the high column density end. We refer to Molinari et al. (2016) for an observed relation between the luminosity-to-mass ratio and the gas temperature of high-mass molecular clumps. The cloud pre-dispersal stage and the cloud dispersed stage may be separated by s₁ ∼ −4 as well. Exactly how the cloud pre-dispersal stage and the cloud dispersed stage are separated from the N-PDF may be sensitive to the cloud initial condition and the subsequent star-formation activities. The more precise loci in the parameter space of the power-law slope of N-PDF, the 2PT function, the N-CCDF, and the luminosity-to-mass ratios, need to be defined with a larger survey. In this sense, the molecular cloud structures of G14.225−0.506 and G28.34+0.06 are on the transition between the cloud integration and the stellar assembly stages, when the massive molecular clumps/cores are forming. The massive clump P2 in G28.34+0.06 already formed high-mass stars, which is a quicker process than the evolution of cloud structures. The wiggling nature of the N-PDF of G14.225−0.506 with a "bump" in between η = 2–3 may be due to the fact that very massive, parsec-scale forming molecular clumps have not yet reached their highest column densities (see also discussion in Section 4.2). W43-south has evolved to a stage where its luminosity-to-mass ratio is higher than W43-main, while its N-PDF is steeper, the N-CCDF is lower, and has less significant embedded massive molecular clump(s), as compared with W43-Main. The OB-cluster-forming region G10.2−0.3 appears to be already dispersed by feedback mechanisms, such that its N-PDF has a steep slope, and N-CCDF is as low as the very young source G11.11−0.12. We note that molecular clouds may have different morphological classes (Lin et al. 2016). The parameters that separate the evolutionary stages may vary with the morphology classes, which needs to be examined with larger samples. However, the present difficulty in assembling a large sample is originated from the distance uncertainty/ambiguity. Based on the VLA survey of NH₃ lines toward 62 high-mass star-forming regions with ${L}_{\mathrm{bol}}\,\sim$ 10⁴ ${L}_{\odot }$, Lu et al. (2014) proposed that the clouds can be separated into the different morphological classes of filaments, concentration, and dispersed and a sub-sample not classified yet. These morphological classes may also be related to the evolutionary stages here proposed, if not purely determined by initial conditions.