KFPA Examinations of Young STellar Object Natal Environments (KEYSTONE): Hierarchical Ammonia Structures in Galactic Giant Molecular Clouds

Table 1 lists the 11 clouds observed by KEYSTONE and their distances. Here, we provide a brief overview of each cloud. For more detailed comparisons between the clouds, see the review by Motte et al. (2018a).

2.1.1. W3

2.1.2. Mon R2

2.1.3. Mon R1 and NGC 2264

2.1.4.  Rosette

2.1.5. M16

2.1.6. M17

2.1.7. W48

2.1.8. Cygnus X

2.1.9. NGC 7538

Data were obtained as part of the KEYSTONE survey, a large project on the GBT that mapped NH₃, HC₅N, HC₇N, HNCO, H₂O, CH₃OH, and CCS emission across 11 GMCs at distances between 0.9 and 3 kpc. Observations were conducted between 2016 October and 2019 March for a total of 356.25 observing hours, including overheads. Table 2 summarizes all observed transitions along with their rest frequencies. The 11 GMCs observed by KEYSTONE were selected from the HOBYS survey. The observing strategy for KEYSTONE targeted all filamentary structures where A_V  > 10 mag in the HOBYS column density maps (see B. Ladjelate et al. 2019, in preparation), which is slightly higher than that used in the GAS survey (A_V  > 7; Friesen et al. 2017). Due to the large amount of foreground and/or background contamination along the line of sight to some of the clouds, this extinction threshold does not have much physical meaning but rather is meant to highlight the densest regions in each cloud. The KEYSTONE observations also exclude parts of M16, M17, and W48 that will be mapped with the GBT by the RAMPS (Hogge et al. 2018).

Observations were made with the GBT's KFPA, which has seven beams arranged in the shape of a hexagon with beam centers separated by ∼95'' on the sky. Following the observational setup used in GAS, each cloud was segmented into 10' × 10' tiles that were observed using on-the-fly mapping and frequency-switching for 11 s of on+off integration time (5.5 s on-source and 5.5 s at reference frequency) per beam (i.e., a total of 77 s when summing over all seven beams) for each resolution element. The 10' × 10' tiles were scanned using on-the-fly mapping, covering the observed region in the R.A. and decl. directions. The row separation (∼13'') and spectrometer dump cadence ensured that each resolution element in the map was sampled by more than three samples in both directions, ensuring Nyquist sampling. Each tile took ∼1.3 hr to complete, with 1–3 tiles observed per session. Table 1 lists the number of tiles completed for each cloud. The survey's completeness, defined as the percentage of the HOBYS maps with A_V  > 10 mag observed by KEYSTONE, ranged from 31% for M17 to 100% for Mon R2 and NGC 7538 (see Table 1).

The telescope's pointing and focus were aligned before mapping each tile to account for changes in the optical performance due to, e.g., temperature- and weather-dependent structural deformations. The KFPA receiver's noise diodes were used to measure the off-source system temperatures for each observing session, which are also temperature and weather dependent. Since each of the KFPA's beams has an independent response (i.e., gain), the Moon was observed at least once per session for flux density calibration, if available. The Moon's large angular size compared to the size of the KFPA beam allowed for beam gains to be calculated from single on-source and off-source observations during each observing session. Figure 1 shows the beam gains for the NH₃ (1,1) spectral windows (IFs 6, 7, and 8) averaged over all observations of the Moon for each polarization. Table 3 displays the final beam gains used for flux density calibration, along with the standard deviation for each average.

Zoom In Zoom Out Reset image size

The GBT's VEGAS backend was configured with eight spectral windows, each 23.44 MHz wide. All five NH₃ transitions (1,1) up to (5,5), along with HC₇N (9–8) and CH₃OH (10₁–9₂) A⁻, were observed in seven of the windows across all seven of the KFPA beams. The eighth VEGAS window covered H₂O (6₁₆–5₂₃), HC₅N (8–7), HC₇N (19–18), HNCO (1_0,1–0_0,0), CH₃OH (12₂–11₁) A⁻, and CCS (2₀–1₀) in only the central KFPA beam. The GBT beam has a FWHM of 32'' at the NH₃ (1,1) rest frequency. This VEGAS configuration is the same as that used by the RAMPS (Hogge et al. 2018).

In this paper, we present the NH₃ (1,1) and (2,2) emission. Other lines will be presented in future KEYSTONE papers. We also identify H₂O (6₁₆–5₂₃) maser emission by eye to include in figures presented in Section 4, but leave the full presentation of those data and a more thorough maser identification technique to G. White et al. (2019, in preparation). The NH₃ data were reduced using gbtpipe, ²⁴ a Python-packaged version of the standard GBT reduction pipeline. The data were calibrated and output as 3D FITS spectral cubes, with R.A., decl., and spectral frequency comprising each axis. The on-the-fly observations were mapped to a grid of square pixels with width of 88, which corresponds to ∼3.5 pixels per FWHM beam of the NH₃ (1,1) line. The spectrum corresponding to each spatial pixel was determined using a weighted average of on-the-fly integrations from all seven beams of the KFPA, including those samples with separations less than one FWHM beam size away from a given map pixel. The weighting scheme was a Gaussian-tapered Bessel function, as described in Friesen et al. (2017) following Mangum et al. (2007). This procedure results in data cubes with a resolution of 32'' and the dense sampling from the mapping strategy and multiple receiver feeds produces high-quality maps without discernible scanning patterns in the image or Fourier domain.

The pipeline also subtracts a first-order polynomial fit to the channels on the edges of each scan prior to gridding to remove any shape in the spectral baselines introduced by, e.g., instrumental effects. The pixel size of the final data cubes is 88, with a spectral resolution of 5.7 kHz, or 0.07 km s⁻¹.

To remove any remaining shape in the spectral baselines, we perform an additional round of per-pixel baseline fitting similar to the method described in Hogge et al. (2018). Namely, a sliding window with a width of 31 channels is used to calculate a "local" standard deviation for every channel in a spectrum. For the 15 channels at each end of the spectrum, the first and last 31 channels are used as the "local" windows, while all other channels are at the center of their "local" window. From the standard deviation distribution of all "local" windows, the central channels belonging to the lowest two quintiles are used for the baseline fit. Thus, channels that belong to an emission line or noise spike are excluded from the baseline fit due to their high "local" standard deviation relative to the non-emission-line channels in the spectrum. Next, polynomials up to third order are fit to the selected channels. A reduced chi-squared value is then calculated for each of the best-fit polynomials against the full spectrum. Finally, the polynomial with the lowest reduced chi-squared value is subtracted from the original spectrum. This baseline subtraction technique is publicly available, ²⁵ along with the full KEYSTONE data reduction code base. The final baseline-subtracted NH₃ (1,1) and (2,2) data cubes are publicly available. ²⁶

The system temperatures for the observations were typically 40–50 K, with a median of ∼43 K. Figure 2 shows histograms of the rms noise for the NH₃ (1,1) and (2,2) maps of each cloud. We calculate the rms using the channels in the best-fit model from our line-fitting procedure (see Section 3.1) with brightness lower than 0.0125 K. While the medians of the rms distributions range from 0.13 to 0.2 K, most of the distributions have a peak below 0.15 K. M17 and M16 have slightly higher noise than the other regions since they are the lowest decl. sources observed (∼−16° and −13°, respectively).

Zoom In Zoom Out Reset image size

Herschel Level 2.5 data products at 70, 160, 250, 350, and 500 μm for each KEYSTONE region were downloaded from the European Space Agency Herschel Science Archive. ²⁷ These maps were originally observed by the HOBYS (Motte et al. 2010) and have spatial resolutions of 84, 135, 182, 249, and 363, respectively. Although the HOBYS team has released dense core and protostar catalogs for Mon R2 (Rayner et al. 2017), W3 (A. Rivera-Ingraham et al. 2019, in preparation), Cygnus X North (Bontemps et al., in preparation), NGC 6334 (Tigé et al. 2017), and NGC 6537 (Russeil et al. 2019), no catalogs have yet been released for many of the clouds targeted by KEYSTONE. In this paper, we use the Herschel 160–500 μm maps to estimate the H₂ column densities and masses of structures identified in the KEYSTONE observations (see Section 3.4). Additionally, we use the 70 μm maps to identify embedded protostars in each cloud (see Section 3.6).

To estimate H₂ column densities for each region, spectral energy distributions (SEDs) were created by combining the 160–500 μm maps for each observed pixel. Full details of the SED-fitting method are described in Singh et al. (2019, in preparation), but are similar to the method applied in all HOBYS papers (see, e.g., B. Ladjelate et al. 2019, in preparation). Here, we provide a brief summary of the process: first, a zero-level offset was added to the 160 μm map based on Planck observations to account for background continuum emission not included in the Herschel data. The Herschel Level 2.5 products for 250–500 μm already have this offset applied, so no additional offsets were added to those maps. Next, all maps were convolved to a resolution of 363 and aligned to the same pixel grid as the 500 μm map. SEDs were then assembled on a pixel-by-pixel basis and a modified blackbody model of the form I_λ  = B_λ (T_D )κ_λ Σ was fit to the data, where I_λ is the surface brightness of the emission, B_λ (T_D ) is the Planck blackbody function at dust temperature T_D , and κ_λ is the dust opacity defined as κ_λ  = 0.1(λ/300 μm)^−β cm² g⁻¹ following Hildebrand (1983) and assuming a gas-to-dust ratio of 100. The dust emissivity, β, varies between 1.2 and 2.0 for each pixel and is based on Planck-derived dust models (Planck Collaboration et al. 2014) that were resampled to the same pixel grid as the Herschel maps. A stacked histogram showing the β distribution across all pixels used for SED fitting in each region is displayed in Figure 3. The β distributions vary from cloud to cloud, with the highest values observed in clouds close to the Galactic plane such as W48, M17, and M16. We used Planck-derived values of β because the Herschel data include only the portion of the dust SED close to the intensity peak. The position of the intensity peak is a function of both T_D and β (e.g., λ_peak = 1.493/T_D (3 + β) for a modified blackbody; Elia & Pezzuto 2016) and it is not possible to remove the degeneracy between these two parameters unless data at longer wavelengths, where I_ν  ∝ λ^−β , are used. Since the Planck data include observations down to 850 μm, they are more capable of constraining β than the Herschel data.

Zoom In Zoom Out Reset image size

The gas surface mass density, Σ, and dust temperature were left as free parameters during the fitting procedure. The resulting best-fit model's Σ was converted to H₂ column density, N(H₂), using Σ = μ_H m_H N(H₂), where μ_H = 2.8 is the mean molecular weight per hydrogen molecule, which assumes the relative mass ratios of hydrogen, helium, and metals are 0.71, 0.27, and 0.02, respectively (see, e.g., Appendix A in Kauffmann et al. 2008), and m_H is the mass of a hydrogen atom. The SED-fitting procedure failed to converge for a small fraction of pixels where the dust continuum emission was saturated. The percentage of affected pixels for the KEYSTONE clouds affected are: M17 (0.09%), W48 (0.008% of pixels), Cygnus X North (0.006% of pixels), NGC 7538 (0.01% of pixels), W3 (0.01% of pixels), and Mon R2 (0.001% of pixels). For the affected pixels, we replace their values with the median column density of the 10 closest pixels with reliable SED fits. As such, this is likely a lower limit to the true column density for those pixels. Any dendrogram-identified leaf (see Section 3.3) that overlaps with one of these affected pixels is also flagged in all catalogs and analyses. The number of leaves in each cloud that include affected pixels are: M17 (1 of 38 leaves), W48 (1 of 100 leaves), Cygnus X North (2 of 200 leaves), NGC 7538 (1 of 73 leaves), and W3 (2 of 84 leaves).

The main difference between the column densities derived in this paper and those of the HOBYS collaboration (B. Ladjelate et al. 2019, in preparation) involves the assumptions on β. Specifically, the HOBYS column density maps assume β = 2 for all pixels while we use 1.2 ≤ β ≤ 2.0 based on Planck dust models that constrain β on large spatial scales (Planck Collaboration et al. 2014). Our lower values of β result in comparatively lower column densities in our maps. For instance, the HOBYS team has released the H₂ column density maps and core/protostar catalog for Mon R2 (Rayner et al. 2017). We find that the Rayner et al. (2017) column densities are on average a factor of ∼2.5 higher than those derived in this paper. Although the higher column densities in the HOBYS maps would lead to larger structure masses in our analysis, we discuss in Section 3.4 that the method used to convert the column densities into structure masses is likely a larger source of uncertainty than the β assumption. Moreover, we also recovered 22 of the 28 (∼79%) protostars identified by Rayner et al. (2017), with the six discrepant sources located in the central Mon R2 hub, which is bright at 70 μm. This suggests that our protostar extraction is likely confusion limited in bright hubs, but can efficiently recover sources that are more isolated.

C¹⁸O (3–2) data cubes observed by the HARP-ACSIS spectrometer on the James Clerk Maxwell Telescope (JCMT) were obtained from the JCMT Science Archive, ²⁸ which is hosted by the Canadian Astronomy Data Centre. Of the 11 clouds observed by KEYSTONE, six were found to have publicly available C¹⁸O (3–2) data cubes in the JCMT Science Archive: Cygnus X North, Cygnus X South, M16, M17, NGC 7538, and W3. The native spectral resolution of the C¹⁸O (3–2) cubes was ∼0.056 km s⁻¹ and the spatial resolution was 153. To match better the spatial and spectral resolution of our NH₃ observations and improve sensitivity, we smoothed the C¹⁸O (3–2) maps to a spatial resolution of 32'' and spectral resolution of 0.11 km s⁻¹. In Section 4.3, we describe how Gaussian line fitting of these data cubes is used to estimate the external, turbulent pressure on the ammonia structures observed by KEYSTONE.

The NH₃ (1,1) and (2,2) lines were used to estimate the excitation temperature (T_ex), kinetic gas temperature (T_K ), centroid velocity (V_LSR), velocity dispersion (σ), and para-NH₃ column density (${N}_{\mathrm{para} \mbox{-} {\mathrm{NH}}_{3}}$) for each pixel. We adopted the line-fitting method of the GAS described in Friesen et al. (2017), which uses the coldammonia model in the pyspeckit Python package (Ginsburg & Mirocha 2011) to generate model ammonia spectra under the assumptions of local thermodynamic equilibrium and a single velocity component along the line of sight. While most of the KEYSTONE spectra are well characterized by a single velocity component, we do see signs of multiple velocity components that are closely separated along the spectral axis in regions of W48 and M17. For those spectra, our single velocity component fitting will produce a best-fit model that has a broadened line width to account for the larger width of the emission line features in the spectrum. In a future KEYSTONE paper, we plan to implement a multiple velocity component fitting method that will robustly identify spectra with more than one velocity component and estimate better the line widths for those spectra (J. Keown et al. 2019, in preparation).

The GAS line-fitting pipeline ²⁹ was applied to all pixels with NH₃ (1,1) signal-to-noise ratio (S/N) > 3, where S/N is measured from the ratio of peak emission-line intensity to the rms of the off-line channels in the spectrum. In addition to the minimum S/N threshold, pixels were excluded from our final parameter maps if they did not meet the following constraints on the best-fit model parameters and uncertainties:

• 1.  
  5 K < T_K  < 40 K (outside this range, the NH₃ (1,1) and (2,2) lines cannot constrain T_K );
• 2.  
  0.05 km s⁻¹ < σ < 2.0 km s⁻¹ (below 0.05 km s⁻¹ is unrealistic since our channel width is only ∼0.07 km s⁻¹; above 2.0 km s⁻¹ is uncharacteristic of NH₃ (1,1) emission in the observed star-forming environments (e.g., Olmi et al. 2010; Pillai et al. 2011) and likely indicates the presence of strong outflows or multiple velocity components along the line of sight);
• 3.  
  ${N}_{\mathrm{para} \mbox{-} {\mathrm{NH}}_{3}}$ < 10¹⁶ cm⁻² (above 10¹⁶ cm⁻² is uncharacteristic of NH₃ emission in the observed star-forming environments e.g., Olmi et al. 2010);
• 4.  
  T_K,err < 5 K;
• 5.  
  σ_err < 2.0 km s⁻¹;
• 6.  
  V_LSR,err < 1 km s⁻¹;
• 7.  
  (log ${N}_{\mathrm{para} \mbox{-} {\mathrm{NH}}_{3}}$)_err < 2 ;

where 4–7 are included to cull fits that were unable to converge.

The final parameter maps for each region are shown in Figures 4–15. To compare each region's ammonia emission to its dust continuum emission, we also plot Herschel H₂ column density contours over the ammonia parameter maps. The ammonia emission tends to occur where total extinction in the V band (A_V ) is larger than ∼6–8 mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comparison of the T_K and σ histograms for each region is presented in Figure 16. Although the T_K distributions are consistent for most of the regions, there are significant temperature differences between the regions with the lowest temperatures (Mon R1 and W3-west) compared to the highest-temperature regions (M17 and Mon R2). Similarly, the σ distributions are fairly consistent across regions, with peak values of 0.3–0.7 km s⁻¹. There are several regions (NGC 7538, W48, M17), however, that have a tail of pixels with large line widths >1 km s⁻¹. These large line-width tails are likely due to a higher fraction of pixels with strong outflows or multiple velocity components along the line of sight.

Zoom In Zoom Out Reset image size

The best-fit models from the NH₃ (1,1) line fitting described in Section 3.1 were used to identify the channels to integrate for producing NH₃ (1,1) integrated intensity maps. Namely, the spectral channels in the best-fit models that were brighter than 0.0125 K were included in the integration. This threshold was selected to include only the channels that are part of an emission line in the best-fit models. Since the off-line channels in a pyspeckit model are slightly above zero due to machine precision, the threshold of 0.0125 K provides a conservative distinction between emission-line and off-line channels in the models. For pixels that did not have any channels above that brightness criterion, we used the set of spectral channels centered on the mean cloud centroid velocity with a range defined by the mean cloud line width. In addition, we blanked all pixels within three pixels from the map edges since they had lower coverage by the KFPA and typically had higher noise. Figures 17–28 show the final NH₃ (1,1) integrated intensity maps for each region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hierarchical nature of molecular clouds (e.g., Falgarone & Puget 1986; Lada 1992; Bonnell et al. 2003) warrants a structure-identification method that handles features with different sizes, shapes, and spatial scales. Dendrograms are a proven identification method that excels at identifying such hierarchical features in both continuum (e.g., Kirk et al. 2013b; Könyves et al. 2015) and molecular line-emission observations (e.g., Rosolowsky et al. 2008; Goodman et al. 2009; Lee et al. 2014; Seo et al. 2015; Friesen et al. 2017; Keown et al. 2017) and simulations (Boyden et al. 2016; Koch et al. 2017; Boyden et al. 2018). This ability arises from the tree-diagram architecture of dendrogram algorithms, which first identifies the pixels in a map that represent local maxima. Next, structures are assembled around the local maxima by joining nearby fainter pixels. These top-level leaves are grown until they either merge with another nearby leaf, at which point they are connected by a branch, or reach a pre-defined noise threshold below which no more pixels are added to the structure. The lowest-level structures above this noise threshold that are connected to branches are known as trunks.

Due to the hyperfine structure of NH₃ (1,1) emission, each hyperfine group would be detected as a distinct structure in a 3D dendrogram extraction of the KEYSTONE spectral cubes. To "remove" the hyperfine structures of NH₃ (1,1), some authors have created a single-Gaussian cube from the line width, peak brightness temperature, and centroid velocity measured from the NH₃ (1,1) emission (e.g., Friesen et al. 2017; Keown et al. 2017). Such a single-Gaussian cube attempts to represent how the ammonia emission would appear without hyperfine splitting. Unless the ammonia emission is fit using a model with multiple velocity components along the line of sight, however, the output single-Gaussian cube does not account for emission with multiple velocity components. Instead, a robust multiple velocity component line-fitting method would first need to be applied to the data to take full advantage of a 3D dendrogram extraction of NH₃ (1,1) cubes. The multiple velocity component models would then allow for the creation of a "multi-Gaussian" cube that removes the hyperfine structures of NH₃ (1,1) while preserving the presence of multiple velocity components along the line of sight. Although a multiple velocity component NH₃ (1,1) line-fitting method has been developed in another KEYSTONE paper (J. Keown et al. 2019, in preparation), the analysis presented here neglects multiple velocity components along the line of sight.

Here, we instead perform a dendrogram analysis of the NH₃ (1,1) integrated intensity maps described in Section 3.2. When the observed emission has only a single velocity component along the line of sight, a 2D dendrogram extraction of the integrated intensity maps will produce similar results as a 3D extraction of the full emission cube. Since the majority of the KEYSTONE observations appear to lack multiple velocity components, a 2D analysis is warranted. We defer a 3D dendrogram analysis of the ammonia data to a future KEYSTONE paper.

The astrodendro Python package was applied to the integrated intensity map for each region. For consistency with the ammonia dendrogram analyses by Friesen et al. (2016) and Keown et al. (2017), we chose the following values for the dendrogram algorithm input parameters.

• 1.  
  min_value = 5 × rms, where rms is the rms noise measured in a region of the integrated intensity map where no emission was detected. For clouds with highly variable noise in the integrated intensity map, the rms was calculated using an emission-free region representative of the highest-noise portion of the map. While this conservative approach may leave some low-brightness sources undetected, it reduces the amount of spurious noise sources detected by the dendrogram. min_value is the lowest intensity a pixel can have to be joined to a neighboring structure.
• 2.  
  min_delta = 2 × rms, where rms is the same as described for min_value. min_delta is the minimum difference in brightness between two structures before they are merged into a single structure.
• 3.  
  min_npix = 10 pixels. min_npix is the minimum number of pixels a structure must contain to remain independent. This parameter prevents noise spikes from being identified as sources.

After running the dendrogram algorithm on the maps, we cull leaf sources from our final catalog that do not meet the following criteria.

• 1.  
  The total area of the leaf, in terms of all the pixels associated with it, must be larger than the total area of the GBT beam. This criterion ensures further that small noise spikes are excluded from our final catalog and analyses.
• 2.  
  The leaf contains at least one pixel that was reliably fit by the NH₃ line-fitting method described in Section 3.1. Here, a reliably fit pixel is one that passes the seven constraints listed in Section 3.1.

Figure 29 shows an example tree diagram for the dendrogram extraction of the Mon R2 region. Leaves that do not pass our selection criteria are shown as black vertical lines, while robust leaves are shown in blue and parent structures are shown in red. The tree diagram shows that our selection criteria preferentially cull isolated leaves that are not associated with larger-scale parent structures and are likely noise spikes in the map. Table 4 provides a sample catalog of the leaves that pass our selection criteria in W3-west. Similar catalogs for all 11 KEYSTONE regions are available online. Of the 970 total leaves identified by the dendrograms in each region, the final catalog includes a total of 856 leaves (∼88%) that passed all of the culling criteria. Figures 17–28 show the final catalog leaf masks overlaid atop the NH₃ (1,1) integrated intensity map for each region. These masks show the full extent of all pixels associated with each leaf. For the remainder of the paper, we refer to leaves and clumps synonymously since the ammonia-identified leaves represent the dense gas structures from which new stars may form.

Zoom In Zoom Out Reset image size

Table 4. W3-west NH₃ (1,1) Leaves Catalog 1

ID	R.A.	Decl.	PA	σmajor	σminor	Reff	Mobs	Mclip	TK	${\sigma }_{{\mathrm{NH}}_{3}}$	${V}_{\mathrm{LSR},{\mathrm{NH}}_{3}}$	log(${N}_{\mathrm{para} \mbox{-} {\mathrm{NH}}_{3}}$)	log(${N}_{{{\rm{H}}}_{2}}$)
	(deg)	(deg)	(deg)	('')	('')	(pc)	(M⊙)	(M⊙)	(K)	(km s−1)	(km s−1)	(cm−2)	(cm−2)
0	34.9844	61.0428	166	16.6	9.9	0.27${}_{-0.06}^{+0.07}$	18.7	3.8	10.9 ± 2.4	0.13 ± 0.02	−15.2	13.7	21.4
1	35.4211	61.0938	167	30.7	12.8	0.41${}_{-0.15}^{+0.27}$	93.7	37.2	13.1 ± 1.7	0.36 ± 0.04	−49.9	13.9	21.5
2	35.3625	61.0890	180	14.6	10.8	0.24${}_{-0.13}^{+0.08}$	23.0	4.2	13.0 ± 2.4	0.36 ± 0.05	−49.7	14.1	21.5
3	35.6144	61.0918	101	14.5	7.0	0.17${}_{-0.12}^{+0.09}$	6.6	1.5	10.2 ± 3.6	0.21 ± 0.04	−49.4	13.6	21.3
4	35.2711	61.1001	196	22.1	12.1	0.33${}_{-0.19}^{+0.22}$	47.5	14.3	12.2 ± 2.5	0.27 ± 0.04	−49.7	14.0	21.5
5	35.5964	61.1042	90	10.8	6.2	0.16${}_{-0.01}^{+0.05}$	8.4	1.3	16.5 ± 3.0	0.34 ± 0.06	−50.1	13.4	21.3
6	35.4752	61.1064	152	17.8	10.6	0.27${}_{-0.07}^{+0.07}$	29.2	7.2	15.4 ± 2.6	0.52 ± 0.07	−50.0	13.8	21.5
7	35.2490	61.1213	63	10.2	8.1	0.17${}_{-0.03}^{+0.02}$	9.7	1.6	13.0 ± 3.9	0.23 ± 0.06	−49.6	13.9	21.5
8	35.2099	61.1515	68	13.8	9.7	0.23${}_{-0.03}^{+0.08}$	14.4	2.3	12.3 ± 3.2	0.25 ± 0.05	−48.8	13.1	21.5
9	35.1733	61.1655	155	14.4	8.2	0.21${}_{-0.03}^{+0.1}$	15.2	4.8	12.8 ± 2.8	0.25 ± 0.05	−48.7	13.8	21.5
10	35.2731	61.4578	86	16.1	9.5	0.22${}_{-0.04}^{+0.13}$	30.8	7.6	14.8 ± 2.5	0.46 ± 0.06	−51.0	13.8	21.5
11	35.2470	61.4522	191	13.4	9.0	0.18${}_{-0.08}^{+0.08}$	16.4	1.2	12.6 ± 3.5	0.67 ± 0.11	−51.4	14.0	21.5

Note. Columns show the following values for each leaf: (1) leaf ID, (2, 3) mean R.A. and decl. in J2000 coordinates, (4) position angle of the major axis, measured in degrees counterclockwise from the west on sky, (5, 6) major and minor axis measured by astrodendro based on the intensity-weighted second moment in the direction of greatest elongation, (7) effective radius defined as R_eff = (A/π)^1/2, where A is the area of all pixels in the leaf's mask on the position–position plane, (8) observed mass of leaf from the sum of its H₂ column density, (9) lower-limit mass of leaf calculated using the "clipping" technique (see the text), (10–13) average kinetic gas temperature, velocity dispersion, NH₃ (1,1) centroid velocity, and para-NH₃ column density for leaf, all weighted by the NH₃ (1,1) integrated intensity map, along with their 1σ uncertainties, (14) median H₂ column density for leaf measured from the spatially filtered column density map and used as N in Equation (7). Both column densities are shown in logarithmic scale. Similar tables for all other KEYSTONE regions are available online. Although the intensity-weighted major and minor axes of some sources are less than the 32'' beam size of the observations, our culling criteria ensure that their total areas when considering all their associated pixels are larger than 32''.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The effective radii of the ammonia-identified leaves were estimated using the area of the leaf masks identified by the dendrogram analysis. Following Kauffmann et al. (2013), we adopt R_eff = (A/π)^1/2 as the effective radius, where A is the area of all pixels in the leaf's mask on the position–position plane. Rosolowsky & Leroy (2006) showed that this area-based radius formulation becomes inaccurate for structures with low S/N and sizes much larger or smaller than the beam size. Although we do enforce that leaves are comprised of pixels with at least 5σ detections (see Section 3.3) and limit structures to being larger than the beam size, R_eff may still be susceptible to such biases. To estimate the uncertainties on our measured radii, we use the method described by Chen et al. (2019a), which uses the radius of the largest circle that fits inside the leaf boundary as the leaf's radius lower limit and the radius of the smallest circle that encompasses the leaf as its radius upper limit. The corresponding uncertainties on R_eff based on these upper and lower limits are listed in Table 4. The uncertainties range from 0.1% to 153% of R_eff, with a median of 35%.

Masses for the ammonia-identified leaves were estimated by summing all the H₂ column density for the pixels inside each leaf's mask. The integrated column densities are then converted to mass assuming the distances to each region listed in Table 1 and a mean molecular weight per hydrogen molecule μ_H = 2.8. In Cygnus X and Mon R2, a small number of leaves (six in Cygnus X North, 14 in Cygnus X South, and one in Mon R2) fall outside the boundaries of our H₂ column density maps. Those sources are, therefore, excluded from our analyses that require a mass determination.

Summing all the column density within the leaf boundaries is likely an upper limit on the mass of the structure. Conversely, a lower limit on the structure's mass can be obtained using the "clipping" technique described in Rosolowsky et al. (2008) and Chen et al. (2019a). Namely, before summing the column density pixels within the leaf boundary, the lowest column density pixel's value is subtracted from all other pixels. This method aims to remove contributions to the structure's observed column density from background sources, but is likely over-estimating the true background contribution in most cases. As such, we adopt the regular integrated column density masses throughout this paper, but show the range the mass could be assuming the "clipped" mass is a lower limit. The clipped masses are also displayed in Table 4 alongside the integrated column density masses. The clipped masses are typically a factor of ∼5 (median) lower than the integrated column density masses.

The left panel of Figure 30 shows the effective radii versus mass for all leaves in our final catalog. A power-law fit to the radius versus mass distribution reveals a best-fit slope of 2.43 ± 0.45, which is consistent with the value of 2 expected for clumps of constant surface density and the value of 3 expected for clumps of constant volume density. The data were fit using a Markov chain Monte Carlo (MCMC) sampler. ³⁰ The MCMC sampling used an orthogonal least-squares likelihood function and uniform priors on the power-law slope and intercept. The best-fit model parameters were taken to be the medians of the accepted parameters in the MCMC chain, while the uncertainty on the best-fit parameters was taken to be the standard deviations of the accepted parameter distributions.

Zoom In Zoom Out Reset image size

To estimate the virial stability of the ammonia-identified leaves, we adopt the virial analysis method described in Keown et al. (2017), which uses the ammonia-derived line widths to derive a virial mass (M_vir) for each structure given by

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}=\displaystyle \frac{5{\sigma }^{2}R}{{aG}}\end{eqnarray} \tag{ 1 }$$

where σ is the velocity dispersion of the core (including both the thermal and nonthermal components), R is the core radius, G is the gravitational constant, and

$$\begin{eqnarray}&&a=\displaystyle \frac{1-k/3}{1-2k/5}\end{eqnarray} \tag{ 2 }$$

is a term which accounts for the radial power-law density profile of a core, where ρ(r) ∝ r^−k (Bertoldi & McKee 1992). M_vir represents the mass that a structure with a given radius and internal kinetic energy would have if it were in virial equilibrium when considering only its gravitational potential and kinetic energies. We also assume that the structure is in a steady state, spherical, isothermal, and has a radial power-law density profile of the form ρ(r) ∝ r^−1.5. Our density profile assumption is motivated by recent observations that found ρ(r) ∝ r^−1.5±0.3 for the inner regions of dense cores (e.g., Pirogov 2009; Kurono et al. 2013) and is likely a more accurate choice than the Gaussian density profile chosen in previous virial analyses (e.g., Keown et al. 2017; Kirk et al. 2017; Pattle et al. 2017). See Keown et al. (2017) for a discussion of the implications of assuming a power-law density profile for sources in a virial analysis. We also set R in Equation (1) to be R_eff and calculate the thermal plus nonthermal velocity dispersion as

$$\begin{eqnarray}&&{\sigma }^{2}={\sigma }_{v}^{2}-\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\mathrm{NH}}_{3}}}+\displaystyle \frac{{k}_{{\rm{B}}}T}{{\mu }_{p}{m}_{{\rm{H}}}}\end{eqnarray} \tag{ 3 }$$

where k_B is Boltzmann's constant, ${m}_{{\mathrm{NH}}_{3}}$ is the molecular mass of NH₃, m_H is the atomic mass of hydrogen, and μ_p is the mean molecular mass of interstellar gas. We use μ_p rather than μ_H for this analysis since μ_p considers the additional contributions of helium, assuming a hydrogen-to-helium abundance ratio of 10 and a negligible admixture of metals, that are required to calculate the thermal gas pressure accurately (2.33; see, e.g., Appendix A in Kauffmann et al. 2008). σ_v and T are the average velocity dispersion and kinetic temperature, respectively, within the core boundaries measured from the NH₃ line-fitting parameter maps. Both the σ_v and T averages are weighted by the NH₃ (1,1) integrated intensity map such that σ_v,avg = w₁ σ₁ + w₂ σ₂⋯ w_n σ_n , where w_n and σ_n are the fraction of the source's integrated intensity and value of the velocity dispersion, respectively, for pixel n.

The ratio of M_vir to the actual observed mass of the structure (M_obs) is known as the virial parameter (α_vir = M_vir/M_obs). This parameter can also be written as ${\alpha }_{\mathrm{vir}}=a2{{\rm{\Omega }}}_{K}/| {{\rm{\Omega }}}_{G}|$, where Ω_K , Ω_G , and a represent the structure's total kinetic energy, gravitational potential energy, and density profile (Equation (2)), respectively (Bertoldi & McKee 1992). Thus, the virial parameter neglects the surface term for the kinetic energy that considers the ambient gas pressure exerted on the structure by the cloud. When α_vir ≥ 2, the structure's internal kinetic energy is large enough to prevent it from being gravitationally bound. Conversely, when α_vir < 2, the structure is deemed gravitationally bound since its gravitational potential energy is large enough relative to its internal kinetic energy (neglecting the effects of magnetic fields and external pressure). Figure 31 shows observed mass versus virial parameter for all leaves in our final catalog. Of the 835 leaves, 523 (∼63%) fall below the α_vir = 2 threshold to be considered gravitationally bound. When looking at each cloud individually, the bound leaf fraction varies from ∼0.3 in Mon R2 to ∼0.9 in M17 and W48. Table 5 lists the virial parameters for the leaves identified in W3-west (similar tables for the other regions are provided online). Table 6 shows the bound leaf fractions for each individual cloud, while Table 7 lists the bound fraction and other population statistics for the full leaf sample.

Zoom In Zoom Out Reset image size

Table 5. W3-west NH₃ (1,1) Leaves Catalog 2

ID	αvir	Mvir	σnt	log $| {{\rm{\Omega }}}_{G}|$	logΩK	log $| {{\rm{\Omega }}}_{\mathrm{Pw}}|$	log $| {{\rm{\Omega }}}_{\mathrm{Pt}}|$	On-filament	hub	Nproto	Bad N(H2) Pixels
		(M⊙)	(km s−1)	(erg)	(erg)	(erg)	(erg)
0	0.67	12.6	0.10	43.8	43.4	44.2	NaN	True	False	0	0
1	0.71	66.3	0.36	45.0	44.7	44.8	NaN	True	False	4	0
2	1.68	38.6	0.36	44.1	44.1	44.1	NaN	True	False	1	0
3	1.82	12.0	0.20	43.1	43.2	43.5	NaN	True	False	0	0
4	0.70	33.3	0.26	44.5	44.2	44.5	NaN	True	False	2	0
5	2.92	24.5	0.33	43.4	43.6	43.5	NaN	True	False	1	0
6	2.72	79.6	0.51	44.2	44.4	44.3	NaN	True	False	0	0
7	1.53	14.9	0.21	43.4	43.4	43.7	NaN	True	False	1	0
8	1.47	21.2	0.24	43.7	43.6	44.0	NaN	True	False	0	0
9	1.32	20.0	0.24	43.8	43.7	43.9	NaN	True	False	1	0
10	1.72	52.9	0.45	44.3	44.4	44.0	NaN	True	False	1	0
11	5.01	82.0	0.67	43.9	44.4	43.8	NaN	True	False	0	0

Note. Columns show the following values for each leaf: (1) leaf ID, (2) virial parameter defined as M_vir/M_obs, (3) virial mass calculated using Equation (1), (4) non-thermal component of the velocity dispersion, (5) gravitational energy density calculated using Equation (5), (6) kinetic energy density calculated using Equation (6), (7) cloud weight pressure energy density calculated using Equation (4), (8) turbulent pressure energy density calculated using Equation (4), with NaN representing a lack of C¹⁸O data for that leaf, (9, 10) whether or not the leaf is on-filament or a hub, (11) number of 70 μm point sources within the leaf's boundary, (12) fraction of pixels in the leaf that were saturated in the H₂ column density map. Similar tables for all other KEYSTONE regions are available online.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We note that the variations in distance to each KEYSTONE target provide a variety of linear scales resolved by our observations. These linear resolution effects are not accounted for with the analysis presented here. In the Appendix, however, we show the impact on the virial parameters of NGC 2264, Mon R1, and Mon R2 (d = 0.9 kpc) if those maps were convolved and downsampled to the linear resolution of the W48 (d = 3.0 kpc) observations. We show that there is indeed a tendency for the identified structures in the distance-adjusted analysis to be bound, which might be affecting the higher fraction of bound structures observed in W48 and M17 (d = 2.0 kpc).

Although dendrograms are able to identify the hierarchical parent structures in which leaves are embedded, they are not optimized for isolating the elongated, filamentary structures that are commonly observed in molecular clouds. To understand how the ammonia-identified leaves in this paper relate to surrounding filamentary structures, we employ a dedicated filament extraction algorithm called getfilaments (Men'shchikov 2013) to identify filaments in each region's H₂ column density map. The getfilaments algorithm is a multi-scale extraction approach designed to identify filamentary background structures in Herschel maps (e.g., Könyves et al. 2015; Marsh et al. 2016; Rivera-Ingraham et al. 2016; Bresnahan et al. 2018). As such, it performs far better than dendrograms at identifying filaments. getfilaments was run on all the Herschel H₂ column density maps using the standard extraction parameters for the algorithm (see Men'shchikov 2013 for the extensive list of getfilaments parameters). The top left panels in Figures 32–43 display the final filament masks, reconstructed up to spatial scales of 145''.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A leaf that has at least one of its pixels corresponding to at least one pixel in a filament mask is designated "on-filament" and all other leaves are termed "off-filament." The "on-filament fraction," i.e., the fraction of leaves in a cloud that are on-filament, is listed in Table 6 and ranges from 0.35 in Cygnus X South to 1.0 in W3-west. For the full sample of 835 leaves with Herschel observations, 454 are on-filament (on-filament fraction of ∼54%).

To compare the number of star-forming leaves in each KEYSTONE region, we identify young, embedded protostars using the Herschel 70 μm maps observed for each cloud. getsources (Men'shchikov et al. 2012), a multi-scale source extraction algorithm designed to identify dense cores and protostars in Herschel observations, was employed to extract point sources at 70 μm only. We adopt the Herschel Gould Belt Survey selection criteria for candidate YSOs described in Section 4.5 of Könyves et al. (2015). The final candidate YSOs are shown as red circles in Figures 32–43. We note that at the distances of the KEYSTONE clouds (0.9 kpc < d < 3.0 kpc), there may be significant incompleteness plus insufficient resolution to separate close sources in the Herschel 70 μm maps. For some regions, there may also be contamination by photodissociation regions that can appear as 70 μm point sources. Since we are only using these 70 μm point sources to indicate which leaves are currently star-forming, rather than using them as a complete catalog of YSOs, the extraction is sufficient for our goals.

We perform a cross-match between the candidate YSO catalogs and leaf catalogs to determine which leaves are protostellar. Leaves with at least one candidate YSO falling within their dendrogram-identified boundary are designated "protostellar." Conversely, leaves without a candidate YSO are termed "starless." The protostellar leaf fraction is listed in Table 6 and ranges from 0.17 in Mon R1 to 0.58 in W3-west. For the 835 leaves with Herschel observations, 288 are protostellar (∼34%).

To understand the impact environment has upon the star formation in the KEYSTONE clouds, we search for relationships between 10 variables of interest: leaf on-filament fraction, "bound" leaf fraction (i.e., the fraction of leaves with α_vir < 2), protostellar leaf fraction, total dense gas mass, total protostar count, total leaf count, dense gas surface mass density, surface protostar density, median cloud kinetic temperature, and cloud distance. The dense gas mass, total protostar count, dense gas surface mass density, and surface protostar density are calculated over the KEYSTONE-mapped boundaries of the cloud where the NH₃ (1,1) integrated intensity is greater than 1.0 K km s⁻¹. The dense gas mass is defined as the integrated H₂ column density within the NH₃ (1,1) 1.0 K km s⁻¹ integrated intensity contour, while the total protostar count is defined as the number of candidate YSOs identified by getsources within that same contour. The threshold of NH₃ (1,1) integrated intensity above 1.0 K km s⁻¹ was chosen since it typically highlights the extent of pixels that were robustly fit during our line-fitting procedure (see, e.g., the lowest NH₃ (1,1) integrated intensity contours in Figures 4–15). Median values of H₂ column density within the NH₃ (1,1) 1.0 K km s⁻¹ integrated intensity contours are typically around 1 × 10²² cm⁻², with a minimum of 4.2 × 10²¹ cm⁻² measured in W3-west and a maximum of 2.3 × 10²² cm⁻² measured in M17.

Figure 44 shows the Pearson correlation coefficients for the matrix of 10 variables. The Pearson correlation coefficient can range from −1 to 1, where −1 and 1 signify that the data fall on a straight line with a negative or positive correlation, respectively. A Pearson correlation coefficient of zero indicates there is no correlation between the variables. For 12 data points and using the Student's t-distribution to test for statistical significance, the null hypothesis that the data have no relationship is rejected at the 99.5% confidence level when the Pearson correlation coefficient is greater than $\sim | 0.71|$. As such, correlations that meet this threshold are outlined by black in Figure 44 and the corresponding scatter plots are shown in Figure 45.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Eight statistically significant correlations were found in the data: (1) decreasing leaf on-filament fraction with increasing dense gas mass, (2) decreasing leaf on-filament fraction with increasing number of protostars, (3) decreasing leaf on-filament fraction with increasing number of leaves, (4) decreasing bound leaf fraction with increasing protostellar surface density, (5) increasing dense gas mass with increasing number of protostars, (6) increasing number of protostars with increasing number of leaves, (7) increasing surface mass density with increasing temperature, and (8) decreasing protostellar surface density with increasing cloud distance.

Since the Pearson correlation coefficient does not take into consideration the uncertainties on each data point, we visualize the scatter of each parameter in Figure 45 by adding errorbars as follows: the lower and upper leaf on-filament fraction errorbars represent the on-filament fractions obtained when using the getfilaments filament masks reconstructed up to spatial scales of 72'' and 290'', respectively, rather than the 145'' scale map. Errorbars for the bound leaf fraction reflect the fractions obtained when assuming the virial parameters for each leaf are at the extremes of their individual uncertainty range. Errorbars for the total protostars, protostellar surface density, and NH₃ leaves represent the $\sqrt{N}$ counting uncertainty. Errorbars for the median kinetic temperature represent the median absolute deviation of each cloud's T_K distribution.

Many of these correlations can be explained with our current understanding of the star formation process. For instance, the positive correlation between dense gas mass and protostars (panel 5 in Figure 45) arises from the well-established relationship between star formation rate (SFR) and dense gas mass (e.g., Gao & Solomon 2004; Wu et al. 2005, 2010; Lada et al. 2012; Stephens et al. 2016; Burkhart 2018). Similarly, the correlation we observe between protostars and total ammonia-identified leaves (Panel 6 in Figure 45) is also due to the SFR–mass relation since the leaves are tracing the dense gas mass in each cloud. Moreover, the negative relationship observed between leaf on-filament fraction and dense gas mass (panel 1 in Figure 45) may also be loosely related to the SFR–mass relation. Since the lower-mass clouds in our sample tend to have higher leaf on-filament fractions, this trend may suggest that the star formation in those environments is more heavily dependent on filaments creating the high densities required to form ammonia leaves. Clouds with higher dense gas mass, however, can form clumps and stars even when filaments are not present due to their more widespread dense gas. The same argument can be applied to the anti-correlations observed between leaf on-filament fraction versus protostars and ammonia-identified leaves (panels 2 and 3 in Figure 45). The SFR–mass relation could also explain the positive relationship between temperature and surface mass density displayed in panel 7 of Figure 45 since a higher SFR could lead to higher gas temperatures. Panel 7 has the lowest Pearson coefficient absolute value (0.71) of all the relationships shown in Figure 45 and is dominated by two data points, however, which suggests it is not as robust as the other relationships presented.

In addition, the negative trend observed between bound leaf fraction and protostellar surface density (panel 4 in Figure 45) may be related to the heating and turbulence injected into the cloud by protostars (e.g., Krumholz & McKee 2008; Hansen et al. 2012; Offner & Chaban 2017; Cunningham et al. 2018; Offner & Liu 2018). As the protostellar density increases, the virial parameters of the leaves may increase due to the higher velocity dispersions and temperatures caused by nearby protostellar (or cluster) feedback (e.g., radiation and outflows). Such a scenario is also suggested by the magnetohydrodynamic simulations of Offner & Chaban (2017), which showed that cores become unbound soon after (<0.1 Myr) the onset of protostar formation due to outflow-induced turbulence. Lastly, the negative correlation between protostar surface density and cloud distance shown in panel 8 is likely related to the larger areas observed for the more distant clouds, which would lower their protostar surface densities. Since protostellar surface density is the only parameter significantly correlated with distance, the distance dependence of the other parameters shown in Figure 44 is likely minimal.

A relationship between dense cores and filamentary structures in dust continuum observations has been noted in several nearby star-forming environments. Polychroni et al. (2013) showed that 67% of the dense cores they identified with CuTEx, a Gaussian fitting and background subtraction source extraction algorithm developed by Molinari et al. (2011), were coincident with filamentary structures in L1641 of Orion A. Using the same filament identification algorithm adopted in this paper, Könyves et al. (2015), Marsh et al. (2016), and J. Di Francesco et al. (2019, in preparation) also found that 75%, 40%, and 40%–80% of starless dense cores in Aquila, Taurus-L1495, and five regions of the Cepheus Flare, respectively, are coincident with filamentary structures. In addition, hydrodynamical simulations of molecular clouds (e.g., Offner et al. 2013) also show a strong correspondence between cores and filamentary structures (Mairs et al. 2014). Although the KEYSTONE clouds analyzed in this paper are located at much greater distances (0.9–3.0 kpc) than L1641 (400 pc), Aquila (250–450 pc), Taurus-L1641 (140 pc), and Cepheus (∼300 pc), we find consistent values for the on-filament fraction (∼0.4–1.0 for KEYSTONE clouds).

Despite the apparent relationship observed between leaves and filaments, we find no significant variations between the virial parameters of the on- and off-filament leaf populations in any of the clouds. For the 454 on-filament leaves, 294 (∼65%) have α_vir < 2, which is consistent with the bound fraction for the entire leaf population (∼63%). Furthermore, Figure 46 shows that the mass, effective radius, average kinetic temperature, and average velocity dispersion for the on-filament and off-filament leaves are essentially identical. Although the more distant clouds (e.g., W48 and NGC 7538) tend to have larger masses and radii than the nearest clouds in our sample (e.g., Mon R1, Mon R2, and NGC 2264; see the Appendix for a discussion of the distance dependence of our results), the similarities between the on-filament and off-filament leaf parameter distributions appear to hold for the individual clouds as well. These similarities indicate that star formation away from filaments might be equally as likely as star formation on filaments in high-mass GMCs since dense gas may be more widespread in those environments. Such a scenario is also suggested by the anti-correlation found between leaf on-filament fraction and dense gas mass discussed in Section 3.7. As dense gas becomes more widespread in high-mass GMCs, the fraction of star formation taking place on filaments may decrease since dense gas is equally as likely to be found away from filaments as it is to be found within them.

Zoom In Zoom Out Reset image size

We also note the existence of a group of ammonia-identified leaves with uncharacteristically larger masses (10²–10³ M_⊙) than the majority of leaves in their respective clouds. Deemed "hubs" or "ridges" based on the nomenclature suggested in Myers (2009), these high-mass leaves are shown by the color green in Figures 32–43. The hubs are located at the intersection of multiple filaments (e.g., Mon R2, NGC 2264, eastern and northern regions in W48) and the ridges are massive filaments (e.g., NGC 7538, M16, southern region of W48). Due to their high masses, these structures all have low virial parameters (α_vir = 0.2–0.5) and are likely gravitationally bound or collapsing. As such, the hubs and ridges may be a result of mass build-up at the locations where filaments are transporting mass from other parts of the cloud.

In addition to the hubs and ridges being coincident with filaments, their mass and size indicate they are likely the precursors of MYSOs and stellar clusters. For example, the right panel of Figure 30 shows the mass–radius plot for the hub/ridge and non-hub/ridge leaves in relation to the Kauffmann & Pillai (2010) empirically derived threshold for massive star formation: m(r) > 870 M_⊙ (r/pc)^1.33. The hubs and ridges tend to be above this threshold, indicating that they will likely form high-mass stars. As shown in Figure 46, the hubs and ridges tend to have higher masses (median log(M_obs) = 3.2 ± 0.5), larger radii (median R_eff = 0.7 ± 0.3), warmer temperatures (median T_K,avg = 19.5 ± 5.1), and larger velocity dispersions (median σ = 0.7 ± 0.1) than the starless leaf population. Instead, the hub and ridge masses, radii, and virial parameters are more similar to the massive star-forming clumps identified by Urquhart et al. (2015), highlighting further their propensity to form massive stars. Furthermore, the hubs and ridges tend to align with the positions of H₂O maser emission (identified by eye and shown as cyan stars in Figures 32–43) also detected by KEYSTONE (G. White et al. 2019, in preparation), which is frequently associated with MYSOs. If dense gas hubs and ridges are indeed the current, or future, sites of MYSO and stellar cluster formation, it would explain the high correspondence observed between those objects and filament intersections (e.g., Myers 2009; Hennemann et al. 2012; Schneider et al. 2012; Li et al. 2016; Motte et al. 2018a).

Figures 32–43 also show the virial parameters of the leaves identified as protostellar and starless in each of the clouds observed by KEYSTONE. In many regions (e.g., Mon R1, Rosette, Cygnus X North, Mon R2, W48), the protostellar core population clearly has larger masses and lower virial parameters than the starless core population. For the 288 protostellar leaves identified, 229 (∼80%) have α_vir < 2. In comparison, 294 of the 547 starless leaves identified (∼54%) have α_vir < 2. As shown in Figure 46, the protostellar population's mass distribution peaks at higher values (median log(mass) = 1.8 ± 0.6 M_⊙) than that of the starless population (median log(mass) = 1.3 ± 0.5 M_⊙), which could explain the lower protostellar virial parameters. In addition, the hubs and ridges identified in the previous section tend to host multiple protostars (see, e.g., NGC 7538, NGC 2264, W48, W3, M16). In several regions, the hubs and ridges host over six protostars. Since the 70 μm maps are typically confusion-limited in the hubs and ridges, their protostar counts are likely under-estimated. This attribute highlights the exceptional environment hubs and ridges provide for cluster formation.

Our final sample of 835 ammonia-identified leaves, for which virial parameters have been measured in a consistent manner, forms one of the largest current samples of dense gas structure virial parameters in high-mass star-forming regions. While many studies have derived virial parameters for single molecular clouds or sub-samples of clumps/cores (e.g., Dunham et al. 2010; Schneider et al. 2010a; Urquhart et al. 2011, 2015; Li et al. 2013; Tan et al. 2013; Friesen et al. 2016; Kirk et al. 2017; Billington et al. 2019), few have viewed virial stability across many clouds. Wienen et al. (2012) observed ammonia emission from a sample of 862 clumps in the inner Galactic disk, but only ∼300 of those had known distance measurements required to estimate masses and virial parameters. Similarly, Kauffmann et al. (2013) compiled a catalog of 1325 virial parameter estimates from previously published catalogs of sources in both high- and low-mass star-forming regions. The Kauffmann et al. catalog, however, featured virial parameters that had been measured using a variety of molecular tracers (¹³CO, NH₃, N₂D⁺), mass estimation methods (dust continuum and near-infrared extinction), and probed varying scales (clouds, clumps, and cores). Nevertheless, we can usefully compare the Kauffmann et al. catalog to our data since it deals primarily with high-mass star-forming regions, uses the same formulation for source effective radius, and uses dust continuum emission to derive the masses for most sources.

Overall, our virial parameters are consistent with those found for the high-mass cores, clumps, and clouds included in the Kauffmann et al. (2013) compilation, which included the clumps observed by Wienen et al. (2012). α_vir ranges from ∼10⁻¹ to 10² and roughly half of the sources fall below the α_vir = 2 threshold for both the Kauffmann et al. (2013) sources and those presented in this paper. This result is in contrast to the virial analyses of Friesen et al. (2016) in Serpens South and Keown et al. (2017) in Cepheus-L1251, which found that nearly all their ammonia-identified leaves had α_vir < 2. Serpens South and Cepheus-L1251 are closer (d ∼ 250–450 pc and 300 pc, respectively) than the clouds observed by KEYSTONE and are thus probing smaller-scale structures, which may explain the higher rate of gravitationally bound leaves in those papers. These results are supported by those of Ohashi et al. (2016) and Chen et al. (2019b), which used ALMA observations of infrared dark clouds to show that α_vir decreases with decreasing spatial scales from filaments to clumps to cores, and may be showing that gravity is more important in the stability of structures at small scales. Alternatively, the KEYSTONE observations may indicate that ammonia is more widespread throughout GMCs than it is in low-mass clouds, producing detectable ammonia emission in both bound and unbound sources. Such a scenario is supported by the observations of Henshaw et al. (2013), which found N₂H⁺ (1–0) to be more extended in the infrared dark cloud G035.39-00.33 than in low-mass star-forming environments.

Observations of more distant cloud clumps (∼2–11 kpc) by the Bolocam Galactic Plane Survey (BGPS; Rosolowsky et al. 2010; Ginsburg et al. 2013), which mapped 1.1 mm dust continuum emission across the Galactic plane, have also indicated low virial parameters (α_vir < 2) for clumps. For instance, Svoboda et al. (2016) combined NH₃ observations with BGPS clump detections to calculate virial parameters for 1640 clumps and found that 76% of starless candidates and 86% of protostellar candidates had α_vir < 2. Similarly, Dunham et al. (2011) found that a separate sample of 456 BGPS clumps had a median α_vir = 0.74. We note, however, that the NH₃ observations used for those analyses were targeted follow-ups to previously identified clumps in the BGPS data. Thus, they may not be tracing the faint NH₃ emission probed by KEYSTONE and included in our leaf catalog. These low-brightness NH₃ sources comprise the lower-mass leaves in our sample that have α_vir > 2 and may be the reason we detect unbound sources that are lacking in the BGPS data.

A similar selection bias for high-brightness sources is likely also impacting clump virial analyses using Herschel Hi-GAL (Molinari et al. 2010) and APEX Telescope ATLASGAL (Schuller et al. 2009) observations. For example, Merello et al. (2019) combined Hi-GAL clump detections with NH₃ catalogs to derive virial parameters for 1068 clumps at distances typically between ∼2 kpc and ∼15 kpc; 72% of the 1068 clumps had 0.1 < α_vir < 1, with a median α_vir of 0.3. A similar virial analysis of 213 Hi-GAL clumps by Traficante et al. (2018) found that 76% have α_vir < 1. Using a sample of 1244 ATLASGAL clump detections with masses much larger (typically >10³ M_⊙) than the leaves presented in this paper, Contreras et al. (2017) found a median α_vir of 1.1. Since these analyses rely on clump catalogs identified by dust continuum observations, however, they likely exclude the faint NH₃ sources detected by KEYSTONE.

Although the virial analysis presented in Section 3.5 compares the gravitational energy of the ammonia-identified leaves with their kinetic energy, it excludes the possible influence of external pressure applied by the leaves' surroundings (Field et al. 2011). For instance, the weight of the molecular cloud in which the leaves are embedded can contribute to their confinement (e.g., Pattle et al. 2015, 2017; Kirk et al. 2017). Here, we add the cloud weight pressure energy density (Ω_Pw) to the virial equation using the technique described in Keown et al. (2017) and Kirk et al. (2017). The three energy densities in the virial equation considered in our analysis are given by the following expressions:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{Pw}}=-4\pi {P}_{w}{R}^{3}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{G}=\displaystyle \frac{-1}{2\sqrt{\pi }}\displaystyle \frac{{{GM}}^{2}}{R}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{K}=\displaystyle \frac{3}{2}M{\sigma }^{2}\end{eqnarray} \tag{ 6 }$$

where M is the observed structure mass, R is the effective radius, G is the gravitational constant, σ² is the same as in Equation (3), and P_w is cloud weight pressure:

$$\begin{eqnarray}&&{P}_{w}=\pi G\bar{N}N{\left({\mu }_{{\rm{H}}}{m}_{{\rm{H}}}\right)}^{2}\end{eqnarray} \tag{ 7 }$$

where $\bar{N}$ is the mean cloud column density and N is the column density at the structure (e.g., McKee 1989; Kirk et al. 2006, 2017). Both $\bar{N}$ and N are measured from a spatially filtered column density map to determine the cloud's mass contribution from large-scale structures. Following other recent virial analyses incorporating turbulent pressure (Kerr et al. 2019; Keown et al. 2017; Kirk et al. 2017), the a trous transform. ³¹ is used to filter spatially each column density map to spatial scales larger than 2ⁿ pixels, where we have chosen n = 4 or 16 pixels for all regions. Figure 47 shows an example spatially filtered column density map for NGC 7538, where 16 pixels is equivalent to ∼96'' or ∼1.3 pc. $\bar{N}$ is calculated as the mean of the spatially filtered column density map, while N represents the mean spatially filtered column density within each leaf's dendrogram-identified boundary.

Zoom In Zoom Out Reset image size

Although we have chosen a single scale for the spatial filtering, a recent virial analysis by Kerr et al. (2019) investigated the impact that varying this scale has upon the cloud weight pressure term. In their analysis of L1688, B18, and NGC 1333, Kerr et al. find that increasing or decreasing the spatial filtering scale by a factor of two results in less than a factor of two difference in the average P_w values for those clouds. As shown below, such a factor is not enough to change our overall conclusions about the virial stability of the observed structures.

The left panel of Figure 48 shows the balance between cloud weight pressure (Ω_Pw), kinetic energy (Ω_K ), and gravitational potential energy (Ω_G ) for the 835 leaves with mass estimates. Leaves to the right of the vertical dotted line are deemed "sub-virial" since their gravitational potential and external pressure are enough to overcome their internal kinetic energy in the absence of magnetic fields. Leaves to the left of the vertical line are "super-virial" since they are not bound by their gravitational and external pressure energy densities. Below the horizontal dotted line are "pressure-dominated" sources that have a higher external pressure energy density than their gravitational energy density. Conversely, "gravity-dominated" sources lie above the horizontal line.

Zoom In Zoom Out Reset image size

As previous virial analyses in nearby low-mass star-forming regions have shown (Pattle et al. 2015, 2017; Keown et al. 2017; Kirk et al. 2017; Chen et al. 2019a; Kerr et al. 2019), most of the leaves are sub-virial and pressure-dominated. In particular, Kerr et al. (2019) showed that 79 of 134 (∼59%) of the cores in their combined sample from the NGC 1333, L1688, and B18 clouds were sub-virial when considering both cloud weight pressure and turbulent pressure. Similarly, we find that ∼69% of the leaves presented in this paper are deemed sub-virial when considering only the cloud weight pressure in the pressure energy density. This implies that the larger H₂ column densities and total cloud masses observed for GMCs are sufficient to create virially bound structures without the necessity of large levels of turbulent pressure that appear to be required in low-mass star-forming environments (see Section 4.5 for a discussion of turbulent pressure).

Although many of the structures are pressure dominated, the large surrounding reservoirs of dense gas in GMCs may facilitate their evolution into the gravitationally dominated regime. For instance, some recent simulations have shown that, even though dense cores may appear as pressure-confined and stable structures at various stages in their evolution, they are still likely to be gaining mass by accreting material from their surroundings and will eventually undergo gravitational collapse (e.g.,"global hierarchical collapse"; Naranjo-Romero et al. 2015; Vázquez-Semadeni et al. 2017, 2019; Ballesteros-Paredes et al. 2018). As such, the observed pressure-dominated state of our ammonia-identified leaves may be a common stage in their evolution from dense gas structures to protostars and clusters. Furthermore, the dearth of structures in the sub-virial and gravitationally dominated regime may indicate that clumps spend the majority of their time being pressure dominated, with a quick transition to being gravitationally dominated and subsequently collapsing to form protostars.

In addition to cloud weight, pressure due to cloud-scale turbulence may have a significant impact on the virial stability of dense cores (e.g., Pattle et al. 2015, 2017; Keown et al. 2017; Kirk et al. 2017; Chen et al. 2019a; Kerr et al. 2019). Here, we calculate the turbulent pressure energy density (Ω_Pt) for a subset of our ammonia-identified leaves following the method described by Keown et al. (2017; see also Pattle et al. 2017 and Kirk et al. 2017 for detailed discussions of turbulent pressure). Ω_Pt is calculated from Equation (4), with P_w being replaced with P_t given by

$$\begin{eqnarray}&&{P}_{T}={\mu }_{{\rm{H}}}{m}_{{\rm{H}}}\times {\rho }_{{{\rm{C}}}^{18}{\rm{O}}}\times {\sigma }_{{{\rm{C}}}^{18}{\rm{O}}}^{2},\end{eqnarray} \tag{ 8 }$$

where ${\sigma }_{{{\rm{C}}}^{18}{\rm{O}}}$ is the velocity dispersion measured from C¹⁸O (3–2), a moderate density tracer, and ${\rho }_{{{\rm{C}}}^{18}{\rm{O}}}$ is the volume density at which the C¹⁸O (3–2) emission originates. Here, we assume the C¹⁸O (3–2) emission is tracing a volume density of 3 × 10⁴ cm⁻³, which is the critical density (n_cr (u−l) = A_ul /γ_ul , where A_ul is the Einstein A coefficient and γ_ul is the collisional rate coefficient for collisions with H₂) of C¹⁸O (3–2) at 20 K (calculated using collisional rate coefficients accessed from the Leiden Atomic and Molecular Database; Schöier et al. 2005).

Six of the 11 KEYSTONE regions were found to have partial JCMT observations of C¹⁸O (3–2) (see Section 2.4 for a discussion of these data and our reduction techniques). A Gaussian model with three free parameters (peak brightness temperature, centroid velocity, and velocity dispersion) was fit to all pixels in the C¹⁸O (3–2) data cubes with S/N > 6 using the nonlinear least-squares curve-fitting method in the scipy.optimize.curvefit Python package. S/N is measured from the ratio of the peak brightness temperature in each spectra to the standard deviation of the off-line spectral channels. The conservative cutoff of S/N > 6 was chosen to remove low-S/N spectra from consideration since they often have higher uncertainties in the fitted parameters. The initial parameter guesses of each fit are based on the brightness and velocity of the peak brightness channel, with a set guess of 1.5 km s⁻¹ for the velocity dispersion. After the line fitting, pixels must meet the following criteria to be included in our final parameter maps:

• 1.  
  σ > 0.05 km s⁻¹ (below 0.05 km s⁻¹ is unrealistic since our channel width is only ∼0.11 km s⁻¹;
• 2.  
  T_peak,err < 1 K;
• 3.  
  σ_err < 0.5 km s⁻¹;
• 4.  
  V_LSR,err < 0.75 km s⁻¹.

The final parameter maps for the velocity dispersion, σ_C18O, are shown in Figures 49 and 50 and . The 52 leaves with at least one reliably fit C¹⁸O (3–2) pixel are shown by red contours in Figures 49 and 50 and . The median value of σ_C18O is calculated within the boundaries of each leaf and converted into a turbulent pressure using Equations (4) and (8). Although most of the C¹⁸O (3–2) spectra are well-characterized by a single Gaussian, some do show wings that may be due to outflows or multiple velocity components along the line of sight. These areas can be seen in the velocity dispersion maps as sharp increases in σ_C18O. Since most of the leaves do not overlap with these sharp transitions, our single Gaussian fit is likely sufficient for our velocity dispersion estimates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right panel of Figure 48 shows the virial plane for structures with C¹⁸O (3–2) measurements when using Ω_Pt as the pressure term in the virial equation. All of these structures fall within the pressure-dominated, sub-virial quadrant. The ratio of Ω_Pt/Ω_Pw for these structures ranges from ∼1 to ∼200, with a median of ∼7. This would suggest that either cloud-scale turbulence or global collapse, rather than cloud weight, is the dominant contributor to the pressure term in the virial equation for these ammonia-identified leaves.

Since the turbulent pressure calculation is sensitive to the assumed value of ${\rho }_{{{\rm{C}}}^{18}{\rm{O}}}$, we also calculate P_T using a factor of 10 lower value (3 × 10³ cm⁻³) for ${\rho }_{{{\rm{C}}}^{18}{\rm{O}}}$ and show the difference as errorbars in Figure 48. This lower density is more characteristic of the effective excitation density of C¹⁸O (3–2), which is often 1–2 orders of magnitude lower than critical densities (Shirley 2015). Under this assumption, the ratio of Ω_Pt/Ω_Pw correspondingly drops by a factor of 10, with a range from ∼0.1 to ∼20 and a median of ∼0.7. In contrast to the scenario using the higher ${\rho }_{{{\rm{C}}}^{18}{\rm{O}}}$ assumption, this new estimation would suggest that cloud weight is dominant over turbulent pressure.

Several recent virial analyses of nearby star-forming regions such as Ophiuchus, B18, and NGC 1333 have shown that turbulent pressure tends to be larger than cloud weight pressure (e.g., Pattle et al. 2015; Kerr et al. 2019). Conversely, cloud weight pressure appears to be larger than turbulent pressure in Orion A (Kirk et al. 2017). We note, however, that these analyses included turbulent pressure measurements across entire clouds. This approach is in contrast to the analysis we present here, which has turbulent pressure measurements for only a small subset of leaves that are generally concentrated on the most active star formation sites in each cloud observed (e.g., DR21 in Cygnus X, W3(OH) and W3 Main in W3, and M17SW in M17) and tend to qualify as hubs or ridges. As such, this biased sample cannot be used to draw generalizations for the full ammonia-identified leaf catalog presented in this paper. Instead, widespread C¹⁸O mapping across the KEYSTONE clouds is required to investigate further the role of turbulent pressure on cloud structure and core dynamics.

In addition to external pressure, the magnetic field may also play an important role in the virial stability of the ammonia-identified leaves. For instance, both observations (e.g., Tan et al. 2013; Pillai et al. 2015) and simulations (e.g., Peters et al. 2011) suggest that magnetic fields are equally as important as turbulence and gravity for high-mass star formation. Although we currently lack large-scale magnetic field measurements for the KEYSTONE clouds, Auddy et al. (2019) have recently developed a "core field structure" model that predicts the magnetic field strength and fluctuation profile using dense gas kinematics. We reserve a detailed analysis of the clump magnetic fields, however, to a future KEYSTONE paper.