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

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

Jared Keown¹, James Di Francesco¹,², Erik Rosolowsky³, Ayushi Singh⁴,⁵, Charles Figura⁶, Helen Kirk¹,², L. D. Anderson⁷, Michael Chun-Yuan Chen¹, Davide Elia⁸, Rachel Friesen⁹, Adam Ginsburg¹⁰, A. Marston¹¹, Stefano Pezzuto⁸, Eugenio Schisano⁸, Sylvain Bontemps¹², Paola Caselli¹³, Hong-Li Liu¹⁴,¹⁵, Steven Longmore¹⁶, Frédérique Motte¹⁷, Philip C. Myers¹⁸, Stella S. R. Offner¹⁹, Patricio Sanhueza²⁰, Nicola Schneider²¹, Ian Stephens¹⁸, and

James Urquhart²²,²³

and

the KEYSTONE collaboration

¹Department of Physics and Astronomy, University of Victoria, Victoria, BC, V8P 5C2, Canada; jkeown@uvic.ca

²NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada

³Department of Physics, University of Alberta, Edmonton, AB, Canada

⁴Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada

⁵Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George St., Toronto, Ontario, M5S 3H8, Canada

⁶Mathematics, Computer Science, and Physics, Wartburg College, Waverly, IA, USA

⁷Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506, USA

⁸INAF-IAPS, via del Fosso del Cavaliere 100, I-00133 Roma, Italy

⁹National Radio Astronomy Observatory, Charlottesville, VA 22903, USA

¹⁰National Radio Astronomy Observatory, 1003 Lopezville Rd, Socorro, NM 87801, USA

¹¹ESAC/ESA, Camino Bajo del Castillo s/n, E-28691 Villanueva de la Cañada, Spain

¹²OASU/LAB-UMR5804, CNRS, Université Bordeaux, F-33615 Pessac, France

¹³Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, D-85748, Garching, Germany

¹⁴Department of Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong

¹⁵Departamento de Astronomía, Universidad de Concepción, Av. Esteban Iturra s/n, Distrito Universitario, 160-C, Chile

¹⁶Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK

¹⁷Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

¹⁸Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

¹⁹Department of Astronomy, The University of Texas, Austin, TX 78712, USA

²⁰National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

²¹I. Physik. Institut, University of Cologne, Zuülpicher Str. 77, D-50937 Cologne, Germany

²²Centre for Astrophysics and Planetary Science, University of Kent, Canterbury, CT2 7NH, UK

²³Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

Received 2019 June 15; revised 2019 August 20; accepted 2019 August 22; published 2019 October 7

## Abstract

We present initial results from the K-band Focal Plane Array Examinations of Young STellar Object Natal Environments survey, a large project on the 100 m Green Bank Telescope mapping ammonia emission across 11 giant molecular clouds at distances of 0.9–3.0 kpc (Cygnus X North, Cygnus X South, M16, M17, Mon R1, Mon R2, NGC 2264, NGC 7538, Rosette, W3, and W48). This data release includes the NH₃ (1,1) and (2,2) maps for each cloud, which are modeled to produce maps of kinetic temperature, centroid velocity, velocity dispersion, and ammonia column density. Median cloud kinetic temperatures range from 11.4 ± 2.2 K in the coldest cloud (Mon R1) to 23.0 ± 6.5 K in the warmest cloud (M17). Using dendrograms on the NH₃ (1,1) integrated intensity maps, we identify 856 dense gas clumps across the 11 clouds. Depending on the cloud observed, 40%–100% of the clumps are aligned spatially with filaments identified in H₂ column density maps derived from spectral energy distribution fitting of dust continuum emission. A virial analysis reveals that 523 of the 835 clumps (~63%) with mass estimates are bound by gravity alone. We find no significant difference between the virial parameter distributions for clumps aligned with the dust-continuum filaments and those unaligned with filaments. In some clouds, however, hubs or ridges of dense gas with unusually high mass and low virial parameters are located within a single filament or at the intersection of multiple filaments. These hubs and ridges tend to host water maser emission, multiple 70 µm detected protostars, and have masses and radii above an empirical threshold for forming massive stars.

**Key words:** ISM: kinematics and dynamics – ISM: structure – stars: formation

**Supporting material:** machine-readable tables

## 1. Introduction

The ubiquity of filaments in star-forming environments was first revealed by continuum observations of nearby (<300 pc), low-mass star-forming molecular clouds, which showed that filaments are present in both quiescent (Miville-Deschênes et al. 2010; Ward-Thompson et al. 2010) and active (André et al. 2010;

Men'schikov et al. (2010) star-forming regions. These results suggest filaments are created during the molecular cloud formation process prior to the onset of star formation, likely as a result of turbulence (Vázquez-Semadeni et al. 2006; Smith et al. 2014a, 2014b; Federrath 2016) and magnetic fields (Hennebelle 2013; Palmeirim et al. 2013; Seifried & Walch 2015).
Furthermore, prestellar cores, the arguably gravitationally bound structures that likely collapse to form stars, are predominantly found along filaments (Könyves et al. 2010, 2015; Marsh et al. 2016). These results provide evidence that the formation and gravitational collapse of filaments is related to the core and star formation processes in low-mass star-forming environments.

Although the study of nearby molecular clouds undoubtedly provides us with a close-up view of the star formation process, such clouds are not representative of the most productive star-forming engines in our Galaxy due to their low abundance of O- and B-type stars and clusters. To observe large samples of high-mass stars (>8 $M_\odot$) and stellar clusters, we must probe giant molecular clouds (GMCs) at distances typically >300 pc from our solar system. While these distant environments require higher spatial resolution and sensitivity, they are more indicative of the majority of clouds in the Galaxy. Similar to nearby clouds, filamentary networks of dense gas are also prevalent throughout GMCs and have been found to be spatially correlated with signposts of high-mass star formation (e.g., Nguyen Luong et al. 2011; Hill et al. 2012b; Motte et al. 2018b). In particular, massive young stellar objects (MYSOs) and embedded stellar clusters appear to be preferentially located at the intersections of multiple filaments seen in dust continuum observations (Myers 2009; Schneider et al. 2010a, 2012; Hennemann et al. 2012; Li et al. 2016; Motte et al. 2018a). The combination of the pervasiveness of filaments throughout molecular clouds with the finding that clusters form at the intersections of multiple filaments motivates the idea that mass flow along filaments provides the localized high-density conditions necessary to form stellar clusters and the MYSOs that form within them (Schneider et al. 2010a; Friesen et al. 2013; Henshaw et al. 2013; Kirk et al. 2013a; Fukui et al. 2015; Motte et al. 2018a).

While dust continuum emission provides a detailed look at the distribution of dense cores and filaments within molecular clouds, it does not provide the gas velocity dispersion measurements required to understand whether or not those structures are gravitationally bound. Rather, observations of dense gas emission from molecules such as $NH_3$ (ammonia) and $N_2H^+$ (diazenylium) are necessary to probe core and filament kinematics. These tracers provide an advantage over commonly observed carbon-based molecules (e.g., CO) for tracing dense gas because they suffer less from freeze-out onto dust grains at the high densities within dense cores (see, e.g., Di Francesco et al. 2007) and they are also typically optically thin with Gaussian-like profiles that allow an easier interpretation of kinematics. In addition, the hyperfine splitting of ammonia emission provides a convenient method for obtaining optical depths. Since the relative heights of the $NH_3$ hyperfine structures are well known in the optically thin limit, optical depths and excitation temperatures can easily be determined by measuring the intensities of the hyperfine components (Ho & Townes 1983). Furthermore, observations of multiple $NH_3$ transitions allow a kinetic gas temperature to be calculated from the relative intensities of the central hyperfine groups in each transition. This line strength relationship serves as a proxy for the distribution of populations within each excited state (Ho et al. 1979), i.e., the kinetic energy over the observed portion of the cloud.

The combination of dense gas kinematics and temperatures with continuum observations provides a way to measure the virial stability of dense cores and filaments (e.g., Friesen et al. 2016;

Keown et al. 2017; Kirk et al. 2017), the dissipation of turbulence from clouds and filaments to cores (“transition to coherence”; Pineda et al. 2010; Chen et al. 2019a), and the flow of gas along or onto filaments (e.g., Schneider et al. 2010a; Friesen et al. 2013; Henshaw et al. 2013; Kirk et al. 2013a). Such measurements can also be used to determine whether dense structures associated with filament intersections are susceptible to gravitational collapse. If so, the structures may be the precursors of future stellar clusters, further linking filament intersections to the star formation process in GMCs.

Recent large surveys have set out to investigate the connection between dense gas kinematics and star formation by observing ammonia emission throughout different regions of the Galaxy. The Green Bank Ammonia Survey (GAS) mapped $NH_3$ emission throughout the nearby Gould Belt molecular clouds ($d < 500$ pc) where $A_\nu > 7$ (e.g., Friesen et al. 2017; Keown et al. 2017; Kirk et al. 2017; Redaelli et al. 2017; Chen et al. 2019a; Kerr et al. 2019). The Galactic plane, which typically excludes nearby ($<3$ kpc) GMCs, has been mapped in ammonia by the Radio Ammonia Mid-Plane Survey (RAMPS; covering $10^\circ < l < 40^\circ$, $-0^\circ.5 < b < +0^\circ.5$; Hogge et al. 2018) and the $H_2O$ Southern Galactic Plane Survey (covering $-70^\circ > l > 30^\circ$, $-0^\circ.5 < b < +0^\circ.5$; Purcell et al. 2012). Similarly, Urquhart et al. (2011, 2015) observed ammonia and water maser emission from ~600 MYSOs and ultra-compact H II regions as part of the Red MSX Source Survey. While these surveys trace the kinematics of the most quiescent and extreme environments in the Galaxy, they do not cover the nearest GMCs producing massive stars.

Here, we present K-band Focal Plane Array (KFPA) Examinations of Young STellar Object Natal Environments (KEYSTONE, PI: J. Di Francesco), a large project on the Green Bank Telescope (GBT) that has mapped $NH_3$ emission in 11 GMCs at intermediate distances ($0.9$ kpc $< d < 3.0$ kpc) using the KFPA receiver and VEGAS spectrometer on the GBT. KEYSTONE targeted GMCs observable from Green Bank that are part of the Herschel OB Young Stars Survey (HOBYS; Motte et al. 2010), which mapped dust continuum emission in all GMCs out to 3 kpc using the *Herschel Space Observatory*. This sample of molecular cloud complexes presented in Motte et al. (2018a; see also Schneider et al. 2011) gives a complete view of high-mass star formation at distances less than 3 kpc. This sample notably contains the Cygnus X molecular complex (Hennemann et al. 2012; Schneider et al. 2016), the M16/M17 complex (Hill et al. 2012b; Tremblin et al. 2013, 2014), the Monoceres complex (Didelon et al. 2015; Rayner et al. 2017), Rosette (Di Francesco et al. 2010; Motte et al. 2010; Schneider et al. 2010b, 2012), W48 (Nguyen Luong et al. 2011; Rygl et al. 2014), the W3/KR140 complex (Rivera-Ingraham et al. 2013, 2015), NGC 7538 (Fallscheer et al. 2013), plus southern regions not presented here (Hill et al. 2012a; Minier et al. 2013; Tigé et al. 2017). Thus, KEYSTONE provides the kinematic counterpart to the HOBYS survey that is required to understand the relationship between dense gas dynamics and massive stars.

This paper, which is the first KEYSTONE publication, provides an initial look at the $NH_3$ (1,1) and (2,2) emission maps observed in each region, catalogs each region's dense gas clumps, estimates the virial stability of those clumps, and compares the spatial distribution of the clumps to the positions of filaments and protostars identified in *Herschel* observations. Dendrograms, tree-diagrams that identify intensity peaks in a
<table><caption>**Table 1**
KEYSTONE Target GMCs</caption><thead><tr><th>Region</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>Distance<br>(kpc)</th><th>Total Mass<br>(M<sub>☉</sub>)</th><th>Total Area<br>(pc<sup>2</sup>)</th><th>Footprints<sup>a</sup><br>Observed</th><th>Completeness<br>A<sub>V</sub> &gt; 10</th></tr></thead><tbody><tr><td>W3</td><td>02:23:22.140</td><td>+61:36:17.432</td><td>2.0 ± 0.1<sup>b</sup></td><td>1.0E5</td><td>2.7E3</td><td>26+</td><td>98%</td></tr><tr><td>Mon R2</td><td>06:08:25.657</td><td>-06:14:32.812</td><td>0.9 ± 0.1<sup>c,d</sup></td><td>4.9E3</td><td>1.4E2</td><td>5+</td><td>100%</td></tr><tr><td>Mon R1</td><td>06:32:32.294</td><td>+10:27:13.335</td><td>0.9 ± 0.1<sup>c,e</sup></td><td>8.8E3</td><td>1.4E2</td><td>5</td><td>98%</td></tr><tr><td>Rosette</td><td>06:33:38.530</td><td>+04:29:10.771</td><td>1.4 ± 0.1<sup>c</sup></td><td>3.2E4</td><td>7.2E2</td><td>15</td><td>85%</td></tr><tr><td>NGC 2264</td><td>06:40:41.339</td><td>+09:25:42.177</td><td>0.9 ± 0.1<sup>c</sup></td><td>1.0E4</td><td>2.2E2</td><td>8+</td><td>99%</td></tr><tr><td>M16</td><td>18:18:38.140</td><td>-13:39:30.050</td><td>1.8 ± 0.5<sup>f</sup></td><td>8.6E4</td><td>5.6E2</td><td>5</td><td>52%</td></tr><tr><td>M17</td><td>18:19:35.479</td><td>-16:19:09.088</td><td>2.0 ± 0.1<sup>g</sup></td><td>5.0E5</td><td>1.9E3</td><td>9</td><td>31%</td></tr><tr><td>W48</td><td>19:00:52.657</td><td>+01:41:55.338</td><td>3.0<sup>h</sup></td><td>1.4E6</td><td>8.8E3</td><td>13+</td><td>60%</td></tr><tr><td>Cygnus X South</td><td>20:33:42.800</td><td>+39:35:41.356</td><td>1.4 ± 0.1<sup>i</sup></td><td>2.2E5</td><td>2.3E3</td><td>43</td><td>84%</td></tr><tr><td>Cygnus X North</td><td>20:37:14.998</td><td>+41:56:04.742</td><td>1.4 ± 0.1<sup>i</sup></td><td>2.7E5</td><td>3.3E3</td><td>36+</td><td>82%</td></tr><tr><td>NGC 7538</td><td>23:14:50.333</td><td>+61:29:04.744</td><td>2.7 ± 0.1<sup>j</sup></td><td>9.3E4</td><td>1.4E3</td><td>17</td><td>100%</td></tr></tbody></table>

Notes. The R.A. and decl. listed are the mid-point of the entire mapped area. The total mass and total area are calculated as the sum of all H₂ column density and area, respectively, mapped in each cloud by Herschel. The completeness represents the percentage of pixels with $A_V > 10$ in the Herschel H₂ column density maps that were observed by KEYSTONE. We assumed an extinction conversion factor of $N_{H_2}/A_V = 0.94 \times 10^{21}$ (Bohlin et al. 1978). The completion percentages for M16, M17, and W48 account for the RAMPS intended coverage of those regions.

<sup>a</sup> Each footprint is $10' \times 10'$. A “+” denotes that a partially completed tile was also observed in that region.

<sup>b</sup> Hachisuka et al. (2006).

<sup>c</sup> Schlafly et al. (2014).

<sup>d</sup> Lombardi et al. (2011).

<sup>e</sup> Baxter et al. (2009).

<sup>f</sup> Bonatto et al. (2006).

<sup>g</sup> Xu et al. (2011).

<sup>h</sup> Rygl et al. (2010).

<sup>i</sup> Rygl et al. (2012).

<sup>j</sup> Moscadelli et al. (2009).

map and determine their hierarchical structure, are used to select dense gas clumps in each cloud. The top-level structures in the dendrogram hierarchy are often called “leaves,” a term that we use synonymously with “clumps” throughout this paper. In Section 2, we describe our GBT observations and data reduction techniques, along with the archival data that were retrieved for our analysis. In Section 3, we outline the methods used to model the NH₃ data, identify NH₃ structures, derive their stability parameters, and compare their spatial distributions to those of dust continuum filaments. In Section 4, we estimate the cloud weight pressure and turbulent pressure exerted on the NH₃ structures. We conclude with a summary of the paper in Section 5 and a discussion of future analyses using the KEYSTONE data in Section 6.

## 2. Observations and Data Reduction

### 2.1. Targets

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

W3 is part of a larger complex located in the Perseus spiral arm that also includes the W4 and W5 molecular clouds (Megeath et al. 2008). The W3 Main, W3(OH), and AFGL 333 regions on the eastern edge of W3 all show signatures of high-mass star formation that may have been triggered by superbubbles from previous generations of star formation (Oey et al. 2005). W3 Main is a particularly popular source for high-mass star formation studies due to its array of H II regions (Tieftrunk et al. 1997; Colley 1980) powered by a cluster of OB stars

(Megeath et al. 1996; Ojha et al. 2004). For instance, Tieftrunk et al. (1998) used NH₃ (1,1) and (2,2) observations of W3 Main and W3(OH) to show that the stellar clusters are littered with cold, dense gas clumps. More recently, Nakano et al. (2017) mapped the AFGL 333 ridge in NH₃ and found evidence for triggered star formation at the edges of the ridge but quiescent (non-triggered) formation in the ridge center. Similarly, Rivera-Ingraham et al. (2011) argued that both triggered and quiescent star formation are required to explain the young stellar object (YSO) population detected in the cloud. More recent large-scale Herschel mapping of W3 by Rivera-Ingraham et al. (2013, 2015) suggested that the triggered star formation was a result of “convergent constructive feedback,” which involves massive stars serving as triggers for subsequent star formation by funneling gas onto a central massive structure.

In this paper, we present the observations of the southwestern half of W3, which includes the small H II region KR 140 (Kallas & Reich 1980), as a separate region that we named W3-west.

#### 2.1.2. Mon R2

Monoceros R2 (Mon R2) is the most distant member of the larger Orion-Monoceros molecular cloud complex, which also includes the Orion A and Orion B clouds. Wilson et al. (2005) contend that Mon R2 and the Orion clouds share a common origin, as evidenced by the alignment of spurs in their CO emission with the Vela supershell. Mon R2 hosts a central reflection nebula with a high stellar volume density (~9000 stars pc⁻³), including several B-type stars (Carpenter et al. 1997). Didelon et al. (2015) estimated that the size of the four main H II regions in Mon R2 range from 0.1 pc for the central ultra-compact H II region, which they suggest is undergoing pressure-driven large-scale collapse, to 0.8 pc for the most
extended classical H II region. Previous NH₃ mapping by Willson & Folch-Pi (1981) and Montalban et al. (1990) have shown that the H II regions are surrounded by dense gas clumps with masses of 1–65 M☉ and kinetic temperatures of 15–30 K. Moreover, recent *Herschel* dust continuum and C¹⁸O observations by Rayner et al. (2017) showed that the gas and dust in Mon R2 has a distinct hub-spoke geometry, with a central hub of protostars and dense cores that may be fed by several connected filaments. The column density probability distribution function from the *Herschel* observations also shows two power-law tails, suggesting both turbulent- and gravity-dominated regimes in Mon R2 (Schneider et al. 2015; Pokhrel et al. 2016).

A *Herschel*-derived sample of 177 dense cores in Mon R2 was published by Rayner et al. (2017). Their masses span 0.084–24 M☉ and their radii span 0.023–0.3 pc. Of the 177 dense cores identified by Rayner et al. (2017), 29 (~16%) were found to be protostellar and 11 had masses >10 M☉.

### 2.1.3. Mon R1 and NGC 2264

The Monoceros OB1 (Mon OB1) GMC includes NGC 2264, one of the most massive star clusters (~1400 members) within 1 kpc of our position in the Galaxy (Dahm 2008; Teixeira et al. 2012; Rapson et al. 2014). Initial CO and CS mapping of the region revealed several outflows associated with the cluster (e.g., Margulis & Lada 1986; Wolf-Chase et al. 1995). Six Herbig-Haro objects have also been detected within this region (Adams et al. 1979; Walsh et al. 1992; Wang et al. 2003). Ammonia mapping by Lang & Willson (1980) and Pagani & Nguyen-Q-Rieu (1987) revealed that the dense gas in NGC 2264 is comprised of two components, each ~0.9 pc in diameter and separated by 0.9 pc, with kinetic temperatures of ~20 K. In addition, Peretto et al. (2006) used more recent observations of dust continuum and molecular line emission to show that several massive clumps in NGC 2264 indicate infall motions and may comprise an intermediate mode of massive star formation.

Just north of NGC 2264 is a more quiescent region of dense gas where a collection of Class 0/I and II objects are forming (Rapson et al. 2014). We henceforth refer to this northern region as “Mon R1,” which it has been referred to in previous literature (Kutner et al. 1979; Ogura 1984). Large-scale CO mapping covering NGC 2264 and Mon R1 by Oliver et al. (1996) revealed that the kinematics of the region are dominated by the Perseus and Local spiral arms.

### 2.1.4. Rosette

The *Rosette* complex is located in the Monoceros constellation south in decl. from Mon OB1, NGC 2264, and Mon R2 (Román-Zúñiga & Lada 2008). The cloud's emission is dominated by NGC 2244, its central OB association of 70 high-mass stars that has created a large H II region (Wang et al. 2008). *Rosette* has been mapped extensively in CO (Blitz & Thaddeus 1980; Blitz & Stark 1986; Schneider et al. 1998; Heyer et al. 2006), which revealed outflows from the massive proto-binary AFGL 961 (Castelaz et al. 1985). Large-scale *Herschel* dust continuum mapping by Di Francesco et al. (2010) revealed 473 dense clumps throughout *Rosette*, 371 being starless and 102 being protostellar, which includes six protostellar massive dense cores and three prestellar massive dense cores with masses between 20 and 40 M☉ (Motte et al. 2010).

Schneider et al. (2010b) also used the *Herschel* observations to show a negative temperature gradient, positive density gradient, and age sequence (more evolved to younger) as distance from the NGC 2244 cluster increases, highlighting the influence of the OB association upon the star formation in the cloud. In addition, Schneider et al. (2012) note that the massive stars and infrared clusters discovered in *Rosette* tend to align with the intersections of dust-identified filaments, providing compelling evidence that massive star formation occurs at the sites of filament mergers.

### 2.1.5. M16

M16, which is also known as the Eagle Nebula, is an H II region located in the Sagittarius spiral arm (Oliveira 2008). The cloud's structure and temperature are influenced by the open cluster NGC 6611 at its center, which contains 52 OB stars (Evans et al. 2005). For example, Hill et al. (2012b) used *Herschel* dust continuum mapping to show there is a clear dust temperature gradient moving away from the NGC 6611 cluster. Tremblin et al. (2014) also showed that the dust-derived column density probability distribution function in M16 has a second peak at high densities, which they attributed to a compressed zone of gas caused by an expanding shell of ionized gas from NGC 6611. In the south of M16 are the famous “Pillars of Creation” or “elephant trunks” imaged with the *Hubble Space Telescope* (Hester et al. 1996) and with *Herschel* (Hill et al. 2012b; Tremblin et al. 2013). The morphology of the Pillars is caused by the ionizing radiation from the central OB stars in M16 (White et al. 1999; Williams et al. 2001; Gritschneder et al. 2010). In addition, recent CO mapping of M16 by Nishimura et al. (2017) revealed a 10 pc diameter cavity of molecular gas near NGC 6611, providing further evidence of the cluster's impact on the star formation in the GMC.

### 2.1.6. M17

M17 (the Omega Nebula) is located south in decl. from M16 by an angular separation of 2°5 (Oliveira 2008). Elmegreen et al. (1979) used CO mapping, however, to show that M17 and M16 form a continuous molecular cloud structure despite their large angular separation, which is a conclusion supported by recent near-infrared imaging (Comerón et al. 2019). Similar to M16, M17 has a central H II region created by an open cluster (NGC 6618) of 53 OB stars (Hoffmeister et al. 2008). While much of the literature is focused on mapping the molecular gas (e.g., Thronson & Lada 1983; Stutzki et al. 1988; Stutzki & Guesten 1990; Pérez-Beaupuits et al. 2015) and dust continuum (e.g., Gatley et al. 1979; Povich et al. 2009) of the M17SW region near NGC 6618, the whole of M17 has recently been mapped in ¹²CO, ¹³CO, and C¹⁸O by Nishimura et al. (2018) and in ¹²CO, ¹³CO, HCO⁺ and HCN by Q. Nguyen Luong et al. (2019, in preparation). M17SW has also been mapped in NH₃ by Lada (1976) and Guesten & Fiebig (1988), which revealed several distinct velocity components in the dense gas and kinetic temperatures of 30–100 K.

### 2.1.7. W48

At 3 kpc (Rygl et al. 2010), W48 is the most distant HOBYS and thus KEYSTONE target. *Herschel* observations of the cloud by Nguyen Luong et al. (2011) revealed numerous H II regions with extended warm dust emission. The IRDC G035.39 −00.33 region in the north of W48 was also found to host 13
high-mass ($M > 20 M_\odot$), compact (diameters of 0.1–0.2 pc), and dense ((2–20) × 10⁵ cm⁻³) cores that could be the precursors of massive stars (Nguyen Luong et al. 2011). Liu et al. (2018) used dust polarization and NH₃ measurements to show that these clumps are likely supported against gravitational collapse by magnetic fields and turbulence. Similarly, Pillai et al. (2011) used interferometric observations of G35.20−1.74 in the east of W48 to show that the cores there were also massive (~9–250 $M_\odot$), dense (>10⁵ cm⁻³), cold (<20 K), and highly deuterated ([NH₂D/NH₃] > 10%), which suggest they are on the verge of forming protoclusters. With several methanol maser emission line detections (Slysh et al. 1995; Minier et al. 2000; Sugiyama et al. 2008; Surcis et al. 2012), which are a signpost of massive stars, it is clear that W48 is an interesting testbed for high-mass star formation studies.

### 2.1.8. Cygnus X

The Cygnus X molecular cloud complex is one of the most active star-forming regions in the nearby Galaxy (Schneider et al. 2016). It hosts over 1800 protostars (Kryukova et al. 2014) and is a favored target for studies of high-mass star formation due to its high concentration of OB associations (e.g., Hanson 2003; Comerón & Pasquali 2012; Wright et al. 2014). The OB associations range in age and size from the young proto-globular cluster Cyg OB2 (Knödléseder 2000; Wright et al. 2014), harboring nearly 100 O-stars, to the slightly older and smaller Cyg OB1, OB3, and OB9 (Uyaniker et al. 2001). It has been mapped extensively in a variety of molecular gas tracers (Wilson & Mauersberger 1990; Schneider et al. 2006, 2010b, 2016; Csengeri et al. 2011a, 2011b; Pillai et al. 2012; Duarte-Cabral et al. 2013, 2014; Dobashi et al. 2014), dust continuum (Motte et al. 2007; Bontemps et al. 2010; Hennemann et al. 2012), and dust polarization (e.g., Ching et al. 2017).

Although previous papers have treated Cygnus X as a single complex (Schneider et al. 2006; Rygl et al. 2012), our observations split the Cygnus X cloud into a north and south region. The choice to treat Cygnus X North and South as separate regions in our analysis is motivated by the observations of Kryukova et al. (2014), which showed that each has distinct luminosity functions and morphological differences indicative of dissimilar star-forming environments. Cygnus X North contains DR21, the massive ridge where a slew of massive stars are forming, including the high-mass core DR21 (OH) (e.g., Mangum et al. 1991; Csengeri et al. 2011a). Previous observations of the DR21 H II region by Guilloteau et al. (1983) mapped the region in NH₃ (1,1), (2,2), (3,3), and (4,4), which revealed absorption in the (1,1) and (2,2) emission that indicates high excitation temperatures ≥100 K. The southern section of Cygnus X is home to DR15, a cluster of ~200 protostars that sits atop a filamentary pillar extended over 10 pc to the south (Rivera-Gálvez et al. 2015).

### 2.1.9. NGC 7538

NGC 7538 is a GMC associated with the Perseus spiral arm (Kun et al. 2008). It harbors several bright H II regions, most notably around the IRS 1–11 sources in its center (Werner et al. 1979; Mallick et al. 2014). Strong outflows have been observed throughout the cloud (Campbell 1984; Scoville et al. 1986; Sandell et al. 2005; Qiu et al. 2011), one of which has

<table><caption>**Table 2**
KEYSTONE Observed Transitions</caption><thead><tr><th>Molecule</th><th>Transition</th><th>Rest Frequency<sup>a</sup><br>(MHz)</th><th>Number of Beams</th></tr></thead><tbody><tr><td>HC<sub>5</sub>N</td><td>(8–7)</td><td>21301.26</td><td>1</td></tr><tr><td>HC<sub>7</sub>N</td><td>(19–18)</td><td>21431.93</td><td>1</td></tr><tr><td>CH<sub>3</sub>OH</td><td>(12<sub>2</sub>–11<sub>1</sub>) A<sup>−</sup></td><td>21550.34</td><td>1</td></tr><tr><td>HNCO</td><td>(1<sub>0</sub><sub>1</sub>–0<sub>0</sub><sub>0</sub>)</td><td>21981.4706(1)</td><td>1</td></tr><tr><td>H<sub>2</sub>O</td><td>(6<sub>16</sub>–5<sub>23</sub>)</td><td>22235.08</td><td>1</td></tr><tr><td>CCS</td><td>(2<sub>0</sub>–1<sub>0</sub>)</td><td>22344.030</td><td>1</td></tr><tr><td>CH<sub>3</sub>OH</td><td>(10<sub>1</sub>–9<sub>2</sub>) A<sup>−</sup></td><td>23444.78</td><td>7</td></tr><tr><td>NH<sub>3</sub></td><td>(1,1)</td><td>23694.4955</td><td>7</td></tr><tr><td>NH<sub>3</sub></td><td>(2,2)</td><td>23722.6336</td><td>7</td></tr><tr><td>NH<sub>3</sub></td><td>(3,3)</td><td>23870.1296</td><td>7</td></tr><tr><td>HC<sub>5</sub>N</td><td>(9–8)</td><td>23963.9010</td><td>7</td></tr><tr><td>NH<sub>3</sub></td><td>(4,4)</td><td>24139.35</td><td>7</td></tr><tr><td>NH<sub>3</sub></td><td>(5,5)</td><td>24532.92</td><td>7</td></tr></tbody></table>

**Note.**

signatures of a massive (~40 $M_\odot$) accreting Class 0 protostar (Sandell et al. 2003). Dust continuum observations covering the IRS 1–11 sources by Reid & Wilson (2005) showed that the bright IR sources are surrounded by massive cold clumps. Herschel observations by Fallscheer et al. (2013) that covered a wider field of view revealed an evacuated ring structure in the east of NGC 7538, with a string of cold clumps detected along the ring's edge. Fallscheer et al. also detected 13 massive ($M > 40 M_\odot$) and cold ($T < 15$ K) clumps that may be starless or contain embedded Class 0 sources, further highlighting the high-mass star-forming potential of the cloud.

Previous ammonia observations in NGC 7538 have been focused primarily on IRS 1, which has shown a slew of rare emission features such as maser emission in H₂O, the nonmetastable <sup>14</sup>NH₃ (10,6), (10,8), (9,8), and (9,6) transitions, and <sup>15</sup>NH₃ (3,3) (Johnston et al. 1989; Hoffman & Seojin Kim 2011; Hoffman 2012), as well as vibrationally excited ammonia (Schilke et al. 1990).

### 2.2. GBT NH₃ Data

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).

<footnote><sup>a</sup> Accessed from Lovas (2004).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Beam gains for the NH₃ (1,1) spectral windows (IFs 6, 7, and 8) averaged over all Moon observations for each feed and polarization.</figcaption></figure>

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.

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

<table><caption>**Table 3**
Beam Gains</caption><thead><tr><th>Beam</th><th>Polarization L</th><th>Polarization R</th></tr></thead><tbody><tr><td>0</td><td>0.979 (0.050)</td><td>0.944 (0.051)</td></tr><tr><td>1</td><td>0.916 (0.137)</td><td>0.868 (0.126)</td></tr><tr><td>2</td><td>0.875 (0.043)</td><td>0.873 (0.044)</td></tr><tr><td>3</td><td>0.785 (0.084)</td><td>0.780 (0.084)</td></tr><tr><td>4</td><td>0.934 (0.066)</td><td>0.805 (0.070)</td></tr><tr><td>5</td><td>0.742 (0.105)</td><td>0.533 (0.074)</td></tr><tr><td>6</td><td>0.876 (0.072)</td><td>0.972 (0.085)</td></tr></tbody></table>

Note. Average beam gains with $1\sigma$ variations shown in parentheses.

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₀,₀), CH₃OH (1₂₂–1₁₁) 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`,<fnref n="24" /> 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 8"/8, 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 8"/8, 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

<footnote n="24">²⁴ https://github.com/GBTSpectroscopy/gbtpipe</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Histograms of the rms noise for the NH₃ (1,1) (left) and NH₃ (2,2) (right) maps of each cloud. The median of the distribution is displayed as a vertical dotted line, with the corresponding value shown in the upper right corner of each panel. The clouds are ordered from top to bottom by increasing median NH₃ (1,1) rms noise.</figcaption></figure>

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,<fnref n="25" /> along with the full KEYSTONE data reduction code base. The final baseline-subtracted NH₃ (1,1) and (2,2) data cubes are publicly available.<fnref n="26" />

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).

## 2.3. Herschel Dust Continuum Data

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.<fnref n="27" /> These maps were originally observed by the HOBYS (Motte et al. 2010) and have spatial resolutions of 8"/4, 13"/5, 18"/2, 24"/9, and 36"/3, 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

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stacked histogram of the dust emissivity, $\beta$, used for SED fitting of each pixel in the Herschel dust continuum maps for all clouds observed by KEYSTONE. The $\beta$ values are from Planck-derived dust models (Planck Collaboration et al. 2014) that have been resampled to the same pixel grid as the Herschel data.</figcaption></figure>

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 36"/3 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

<footnote n="25">²⁵ https://github.com/GBTAmmoniaSurvey/keystone</footnote>

<footnote n="26">²⁶ doi:10.11570/19.0074.</footnote>

<footnote n="27">²⁷ http://archives.esac.esa.int/hsa/whsa/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Kinetic temperature (top) and NH₃ (1,1) velocity dispersion (bottom) derived from NH₃ (1,1) and (2,2) line fitting of the W3 observations. The black contours show the NH₃ (1,1) integrated intensity at 1.0, 3.5, and 10 K km s⁻¹. The solid and dotted gray contours outline H₂ column densities of 2.8 × 10²¹ cm⁻² and 9.4 × 10²¹ cm⁻², respectively, which are equivalent to a total extinction in the V band of A_V = 3 mag and 10 mag. The 32" beam size is shown as a black dot in the upper left corner of each plot. The thick black and gray lines outline the KEYSTONE and Herschel mapping boundaries, respectively.</figcaption></figure>

the form $I_{\lambda} = B_{\lambda}(T_{D})\kappa_{\lambda} \Sigma$ was fit to the data, where $I_{\lambda}$ is the surface brightness of the emission, $B_{\lambda}(T_{D})$ is the Planck blackbody function at dust temperature $T_{D}$, and $\kappa_{\lambda}$ is the dust opacity defined as $\kappa_{\lambda} = 0.1(\lambda/300~\mu\text{m})^{-\beta} \text{cm}^{-2}\text{g}^{-1}$ following Hildebrand (1983) and assuming a gas-to-dust ratio of 100. The dust emissivity, $\beta$, 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 $\beta$ distribution across all pixels used for SED fitting in each region is displayed in Figure 3. The $\beta$ 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 $\beta$ 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 $\beta$ (e.g., $\lambda_{\text{peak}} = 1.493/T_{D}(3 + \beta)$ for

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Figure 4 for W3-west.</figcaption></figure>

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_{\nu} \propto \lambda^{-\beta}$, are used. Since the Planck data include observations down to 850 µm, they are more capable of constraining $\beta$ than the Herschel data.

The gas surface mass density, $\Sigma$, and dust temperature were left as free parameters during the fitting procedure. The resulting best-fit model's $\Sigma$ was converted to H₂ column density, $N(H_2)$, using $\Sigma = \mu_H m_H N(H_2)$, where $\mu_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
<figure><img src="image_8.png" /><figcaption>Figure 6. Same as Figure 4 for Mon R2.</figcaption></figure>

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 $\beta$. Specifically, the HOBYS column density maps assume $\beta = 2$ for all pixels while we use $1.2 \le \beta \le 2.0$ based on Planck dust models that constrain $\beta$ on large spatial scales (Planck Collaboration et al. 2014). Our lower values of $\beta$ result in comparatively lower column densities in our maps. For instance, the HOBYS team has released the H$_2$ 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 $\beta$ 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 $\mu$m. This suggests that our protostar extraction is likely confused limited in bright hubs, but can efficiently recover sources that are more isolated.

## 2.4. JCMT C$^{18}$O Data

C$^{18}$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,<fnref n="28" /> which is hosted by the Canadian Astronomy Data Centre. Of the 11 clouds observed by KEYSTONE, six were found to have publicly available C$^{18}$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$^{18}$O (3–2) cubes was ~0.056 km s$^{-1}$ and the spatial resolution was 15"/3. To match better the spatial and spectral resolution of our NH$_3$ observations and improve sensitivity, we smoothed the C$^{18}$O (3–2) maps to a

<footnote n="28">²⁸ http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7. Same as Figure 4 for Mon R1.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Same as Figure 4 for Rosette.</figcaption></figure>

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.

# 3. Analysis and Results

## 3.1. NH₃ Line Fitting

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 ($\sigma$), and para-NH₃ column density ($N_{para-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<fnref n="29" /> 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\,\text{K} < T_K < 40\,\text{K}$ (outside this range, the NH₃ (1,1) and (2,2) lines cannot constrain $T_K$);

2. $0.05\,\text{km s}^{-1} < \sigma < 2.0\,\text{km s}^{-1}$ (below $0.05\,\text{km s}^{-1}$ is unrealistic since our channel width is only ~0.07 km s⁻¹; above $2.0\,\text{km s}^{-1}$ 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

<footnote n="29">²⁹ Available at http://gas.readthedocs.io/.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 9. Same as Figure 4 for NGC 2264.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Same as Figure 4 for M16. The solid and dotted gray contours outline H₂ column densities of 5.6 × 10²¹ cm⁻² and 9.4 × 10²¹ cm⁻², respectively, which are equivalent to a total extinction in the V band of A_V = 6 mag and 10 mag.</figcaption></figure>

strong outflows or multiple velocity components along the line of sight);

3. $N_{\text{para-NH}_3} < 10^{16} \text{ cm}^{-2}$ (above $10^{16} \text{ cm}^{-2}$ is uncharacteristic of NH}_3$ emission in the observed star-forming environments e.g., Olmi et al. 2010);

4. $T_{K,\text{err}} < 5 \text{ K}$;

5. $\sigma_{\text{err}} < 2.0 \text{ km s}^{-1}$;

6. $V_{\text{LSR, err}} < 1 \text{ km s}^{-1}$;

7. ($\log N_{\text{para-NH}_3}_{\text{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.

A comparison of the $T_K$ and $\sigma$ 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
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Same as Figure 4 for M17. The solid and dotted gray contours outline H₂ column densities of 7.5 × 10²¹ cm⁻² and 9.4 × 10²¹ cm⁻², respectively, which are equivalent to a total extinction in the V band of A_V = 8 mag and 10 mag.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as Figure 10 for W48.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Same as Figure 4 for Cygnus X South. The solid and dotted gray contours outline H₂ column densities of 4.7 × 10²¹ cm⁻² and 9.4 × 10²¹ cm⁻², respectively, which are equivalent to a total extinction in the V band of A_V = 5 mag and 10 mag.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Same as Figure 13 for Cygnus X North.</figcaption></figure>

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.

## 3.2. NH₃ (1,1) Integrated Intensity Maps

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.
<figure><img src="image_17.png" /><figcaption>Figure 15. Same as Figure 4 for NGC 7538.</figcaption></figure>

Figures 17–28 show the final NH₃ (1,1) integrated intensity maps for each region.

### 3.3. Identifying NH₃ Structures with Dendrograms

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.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Histograms of kinetic temperature (left) and NH₃ (1,1) velocity dispersion (right) for the reliably fit pixels in each region. The median and median absolute deviation of each distribution is printed in the top right corner of each panel. The clouds are ordered from top to bottom by increasing median kinetic temperature.</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>**Figure 17.** NH₃ (1,1) integrated intensity map for W3. Green contours outline leaves identified by a dendrogram analysis of the map that passed the culling criteria listed in Section 3.3.</figcaption></figure>

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.

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Same as Figure 17 for W3-west.</figcaption></figure>

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
<figure><img src="image_22.png" /><figcaption>Figure 19. Same as Figure 17 for Mon R2.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 20. Same as Figure 17 for Mon R1.</figcaption></figure>

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.

### 3.4. Determining Leaf Radii and Masses

The effective radii of the ammonia-identified leaves were estimated using the area of the leaf masks identified by the

<figure><img src="image_24.png" /><figcaption>Figure 21. Same as Figure 17 for Rosette.</figcaption></figure>

dendrogram analysis. Following Kauffmann et al. (2013), we adopt $R_{\text{eff}} = (A/\pi)^{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_{\text{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_{\text{eff}}$ based on these upper and lower limits are listed in Table 4. The uncertainties range from 0.1% to 153% of $R_{\text{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
<table><caption>**Table 4**
W3-west NH₃ (1,1) Leaves Catalog 1</caption><thead><tr><th>ID</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>PA<br>(deg)</th><th>σ<sub>major</sub><br>(")</th><th>σ<sub>minor</sub><br>(")</th><th>R<sub>eff</sub><br>(pc)</th><th>M<sub>obs</sub><br>(M<sub>☉</sub>)</th><th>M<sub>clip</sub><br>(M<sub>☉</sub>)</th><th>T<sub>K</sub><br>(K)</th><th>σ<sub>NH<sub>3</sub></sub><br>(km s<sup>-1</sup>)</th><th>V<sub>LSR,NH<sub>3</sub></sub><br>(km s<sup>-1</sup>)</th><th>log(N<sub>para-NH<sub>3</sub></sub>)<br>(cm<sup>-2</sup>)</th><th>log(N<sub>H<sub>2</sub></sub>)<br>(cm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>0</td><td>34.9844</td><td>61.0428</td><td>166</td><td>16.6</td><td>9.9</td><td>0.27<sup>+0.07</sup><sub>-0.06</sub></td><td>18.7</td><td>3.8</td><td>10.9 ± 2.4</td><td>0.13 ± 0.02</td><td>-15.2</td><td>13.7</td><td>21.4</td></tr><tr><td>1</td><td>35.4211</td><td>61.0938</td><td>167</td><td>30.7</td><td>12.8</td><td>0.41<sup>+0.27</sup><sub>-0.15</sub></td><td>93.7</td><td>37.2</td><td>13.1 ± 1.7</td><td>0.36 ± 0.04</td><td>-49.9</td><td>13.9</td><td>21.5</td></tr><tr><td>2</td><td>35.3625</td><td>61.0890</td><td>180</td><td>14.6</td><td>10.8</td><td>0.24<sup>+0.08</sup><sub>-0.13</sub></td><td>23.0</td><td>4.2</td><td>13.0 ± 2.4</td><td>0.36 ± 0.05</td><td>-49.7</td><td>14.1</td><td>21.5</td></tr><tr><td>3</td><td>35.6144</td><td>61.0918</td><td>101</td><td>14.5</td><td>7.0</td><td>0.17<sup>+0.09</sup><sub>-0.12</sub></td><td>6.6</td><td>1.5</td><td>10.2 ± 3.6</td><td>0.21 ± 0.04</td><td>-49.4</td><td>13.6</td><td>21.3</td></tr><tr><td>4</td><td>35.2711</td><td>61.1001</td><td>196</td><td>22.1</td><td>12.1</td><td>0.33<sup>+0.22</sup><sub>-0.19</sub></td><td>47.5</td><td>14.3</td><td>12.2 ± 2.5</td><td>0.27 ± 0.04</td><td>-49.7</td><td>14.0</td><td>21.5</td></tr><tr><td>5</td><td>35.5964</td><td>61.1042</td><td>90</td><td>10.8</td><td>6.2</td><td>0.16<sup>+0.05</sup><sub>-0.01</sub></td><td>8.4</td><td>1.3</td><td>16.5 ± 3.0</td><td>0.34 ± 0.06</td><td>-50.1</td><td>13.4</td><td>21.3</td></tr><tr><td>6</td><td>35.4752</td><td>61.1064</td><td>152</td><td>17.8</td><td>10.6</td><td>0.27<sup>+0.07</sup><sub>-0.07</sub></td><td>29.2</td><td>7.2</td><td>15.4 ± 2.6</td><td>0.52 ± 0.07</td><td>-50.0</td><td>13.8</td><td>21.5</td></tr><tr><td>7</td><td>35.2490</td><td>61.1213</td><td>63</td><td>10.2</td><td>8.1</td><td>0.17<sup>+0.02</sup><sub>-0.03</sub></td><td>9.7</td><td>1.6</td><td>13.0 ± 3.9</td><td>0.23 ± 0.06</td><td>-49.6</td><td>13.9</td><td>21.5</td></tr><tr><td>8</td><td>35.2099</td><td>61.1515</td><td>68</td><td>13.8</td><td>9.7</td><td>0.23<sup>+0.08</sup><sub>-0.03</sub></td><td>14.4</td><td>2.3</td><td>12.3 ± 3.2</td><td>0.25 ± 0.05</td><td>-48.8</td><td>13.1</td><td>21.5</td></tr><tr><td>9</td><td>35.1733</td><td>61.1655</td><td>155</td><td>14.4</td><td>8.2</td><td>0.21<sup>+0.1</sup><sub>-0.03</sub></td><td>15.2</td><td>4.8</td><td>12.8 ± 2.8</td><td>0.25 ± 0.05</td><td>-48.7</td><td>13.8</td><td>21.5</td></tr><tr><td>10</td><td>35.2731</td><td>61.4578</td><td>86</td><td>16.1</td><td>9.5</td><td>0.22<sup>+0.13</sup><sub>-0.04</sub></td><td>30.8</td><td>7.6</td><td>14.8 ± 2.5</td><td>0.46 ± 0.06</td><td>-51.0</td><td>13.8</td><td>21.5</td></tr><tr><td>11</td><td>35.2470</td><td>61.4522</td><td>191</td><td>13.4</td><td>9.0</td><td>0.18<sup>+0.08</sup><sub>-0.08</sub></td><td>16.4</td><td>1.2</td><td>12.6 ± 3.5</td><td>0.67 ± 0.11</td><td>-51.4</td><td>14.0</td><td>21.5</td></tr></tbody></table>

**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/\pi)^{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\sigma$ 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''.

(This table is available in its entirety in machine-readable form.)

<figure><img src="image_25.png" /><figcaption>**Figure 22.** Same as Figure 17 for NGC 2264.</figcaption></figure>

$\mu_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.<fnref n="30" /> 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.

<footnote n="30">³⁰ The emcee package: https://emcee.readthedocs.io.</footnote>
<figure><img src="image_26.png" /><figcaption>Figure 23. Same as Figure 17 for M16.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 24. Same as Figure 17 for M17.</figcaption></figure>

### 3.5. Virial Analysis

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

$$M_{vir} = \frac{5\sigma^2R}{aG} \qquad (1)$$

where $\sigma$ 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

$$a = \frac{1 - k/3}{1 - 2k/5} \qquad (2)$$

is a term which accounts for the radial power-law density profile of a core, where $\rho(r) \propto 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 $\rho(r) \propto r^{-1.5}$. Our density profile assumption is motivated by recent observations that found $\rho(r) \propto r^{-1.5 \pm 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

$$\sigma^2 = \sigma_v^2 - \frac{k_B T}{m_{NH_3}} + \frac{k_B T}{\mu_p m_H} \qquad (3)$$

where $k_B$ is Boltzmann's constant, $m_{NH_3}$ is the molecular mass of NH$_3$, $m_H$ is the atomic mass of hydrogen, and $\mu_p$ is the mean molecular mass of interstellar gas. We use $\mu_p$ rather than $\mu_H$ for this analysis since $\mu_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). $\sigma_v$ and T are the average velocity dispersion and kinetic temperature, respectively, within the core boundaries measured from the NH$_3$ line-fitting parameter
<figure><img src="image_28.png" /><figcaption>Figure 25. Same as Figure 17 for W48.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 26. Same as Figure 17 for Cygnus X South.</figcaption></figure>

maps. Both the $\sigma_v$ and T averages are weighted by the NH$_3$ (1,1) integrated intensity map such that $\sigma_{v,avg} = w_1\sigma_1 + w_2\sigma_2 + \dots + w_n\sigma_n$, where $w_n$ and $\sigma_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 ($\alpha_{vir} = M_{vir}/M_{obs}$). This parameter can also be written as $\alpha_{vir} = a2\Omega_K/|\Omega_G|$, where $\Omega_K$, $\Omega_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 $\alpha_{vir} \ge 2$, the structure's internal kinetic energy is large enough to prevent it from being gravitationally bound. Conversely, when $\alpha_{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 $\alpha_{vir} = 2$ threshold to be

<figure><img src="image_30.png" /><figcaption>Figure 27. Same as Figure 17 for Cygnus X North.</figcaption></figure>

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.

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
<figure><img src="image_31.png" /><figcaption>**Figure 28.** Same as Figure 17 for NGC 7538.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 29.** Dendrogram tree diagram for Mon R2 showing the peak intensity for each structure identified. Leaves are shown in blue and branches are shown in red. Leaves that were culled based on the selection criteria described in Section 3.3 are shown in black. Leaves that are red denote culled leaves that are part of a branch with an accepted leaf.</figcaption></figure>

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).

### 3.6. Identifying Filaments and Candidate YSOs

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'schikov 2013) to identify filaments in each region's $H_2$ 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_2$ column density maps using the standard extraction parameters for the algorithm (see Men'schikov 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".
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Left: effective radius vs. mass for the leaves identified in each region. The solid black line shows the best power-law fit to the data using an MCMC sampler. The gray lines show 1000 random selections from the full MCMC chain. The dotted line shows a power-law slope of 1.9 (Larson 1981) for comparison. Right: effective radius vs. mass for hubs (green) and non-hubs (black) identified in each region (see Section 4.1 for a discussion of hubs). The dashed black line denotes the empirically derived threshold for massive star formation determined by Kauffmann & Pillai (2010). The red dots with errorbars show the median mass and radius uncertainties in different bins along each axis. All other errorbars show leaves that have mass lower limits due to saturated pixels in the H₂ column density map (see Section 2.3).</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 31.** Virial parameter ($\alpha_{vir}$) vs. mass for the leaves identified in each region. The dashed line denotes $\alpha_{vir} = 2$ when assuming a power-law density profile for the structures. Above this line, structures are deemed to be gravitationally unbound in the absence of magnetic fields or external pressure. The red dots with errorbars show the median mass and radius uncertainties in different bins along each axis. All other errorbars show leaves that have mass lower limits due to saturated pixels in the H₂ column density map (see Section 2.3).</figcaption></figure>

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'schikov 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
<figure><img src="image_35.png" /><figcaption>Figure 33. Same as Figure 32 for W3-west.</figcaption></figure>

<table><caption>Table 5
W3-west NH₃ (1,1) Leaves Catalog 2</caption><thead><tr><th>ID</th><th>α<sub>vir</sub></th><th>M<sub>vir</sub><br>(M<sub>☉</sub>)</th><th>σ<sub>nt</sub><br>(km s<sup>-1</sup>)</th><th>log |Ω<sub>G</sub>|<br>(erg)</th><th>logΩ<sub>K</sub><br>(erg)</th><th>log |Ω<sub>Pw</sub>|<br>(erg)</th><th>log |Ω<sub>Pt</sub>|<br>(erg)</th><th>On-filament</th><th>hub</th><th>N<sub>proto</sub></th><th>Bad N(H<sub>2</sub>) Pixels</th></tr></thead><tbody><tr><td>0</td><td>0.67</td><td>12.6</td><td>0.10</td><td>43.8</td><td>43.4</td><td>44.2</td><td>NaN</td><td>True</td><td>False</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0.71</td><td>66.3</td><td>0.36</td><td>45.0</td><td>44.7</td><td>44.8</td><td>NaN</td><td>True</td><td>False</td><td>4</td><td>0</td></tr><tr><td>2</td><td>1.68</td><td>38.6</td><td>0.36</td><td>44.1</td><td>44.1</td><td>44.1</td><td>NaN</td><td>True</td><td>False</td><td>1</td><td>0</td></tr><tr><td>3</td><td>1.82</td><td>12.0</td><td>0.20</td><td>43.1</td><td>43.2</td><td>43.5</td><td>NaN</td><td>True</td><td>False</td><td>0</td><td>0</td></tr><tr><td>4</td><td>0.70</td><td>33.3</td><td>0.26</td><td>44.5</td><td>44.2</td><td>44.5</td><td>NaN</td><td>True</td><td>False</td><td>2</td><td>0</td></tr><tr><td>5</td><td>2.92</td><td>24.5</td><td>0.33</td><td>43.4</td><td>43.6</td><td>43.5</td><td>NaN</td><td>True</td><td>False</td><td>1</td><td>0</td></tr><tr><td>6</td><td>2.72</td><td>79.6</td><td>0.51</td><td>44.2</td><td>44.4</td><td>44.3</td><td>NaN</td><td>True</td><td>False</td><td>0</td><td>0</td></tr><tr><td>7</td><td>1.53</td><td>14.9</td><td>0.21</td><td>43.4</td><td>43.4</td><td>43.7</td><td>NaN</td><td>True</td><td>False</td><td>1</td><td>0</td></tr><tr><td>8</td><td>1.47</td><td>21.2</td><td>0.24</td><td>43.7</td><td>43.6</td><td>44.0</td><td>NaN</td><td>True</td><td>False</td><td>0</td><td>0</td></tr><tr><td>9</td><td>1.32</td><td>20.0</td><td>0.24</td><td>43.8</td><td>43.7</td><td>43.9</td><td>NaN</td><td>True</td><td>False</td><td>1</td><td>0</td></tr><tr><td>10</td><td>1.72</td><td>52.9</td><td>0.45</td><td>44.3</td><td>44.4</td><td>44.0</td><td>NaN</td><td>True</td><td>False</td><td>1</td><td>0</td></tr><tr><td>11</td><td>5.01</td><td>82.0</td><td>0.67</td><td>43.9</td><td>44.4</td><td>43.8</td><td>NaN</td><td>True</td><td>False</td><td>0</td><td>0</td></tr></tbody></table>

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<sup>18</sup>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<sub>2</sub> column density map. Similar tables for all other KEYSTONE regions are available online.

(This table is available in its entirety in machine-readable form.)
<figure><img src="image_36.png" /><figcaption>**Figure 34.** Same as Figure 32 for Mon R2. Magenta contours on the left side of the upper left panel denote leaves that were not included in our virial analysis since their masses could not be estimated due to their locations being outside our H₂ column density map boundaries.</figcaption></figure>

<table><caption>**Table 6**
Cloud Statistics</caption><thead><tr><th>Cloud</th><th><i>T</i><sub>med</sub><br>(K)</th><th><i>σ</i><sub>med</sub><br>(km s<sup>-1</sup>)</th><th>Protostars</th><th>Area<br>(pc<sup>2</sup>)</th><th>Mass<br>(M<sub>☉</sub>)</th><th>Total Leaves</th><th>Filament Leaves<br>Fraction</th><th>Proto Leaves<br>Fraction</th><th>Bound Leaves<br>Fraction</th></tr></thead><tbody><tr><td>W3</td><td>16.0 ± 4.6</td><td>0.5 ± 0.3</td><td>80</td><td>26.5</td><td>5999</td><td>84</td><td>0.62</td><td>0.38</td><td>0.48</td></tr><tr><td>W3-west</td><td>13.4 ± 2.8</td><td>0.3 ± 0.1</td><td>10</td><td>2.9</td><td>303</td><td>12</td><td>1.00</td><td>0.58</td><td>0.75</td></tr><tr><td>Mon R2</td><td>17.7 ± 4.9</td><td>0.3 ± 0.2</td><td>23</td><td>2.7</td><td>715</td><td>41</td><td>0.73</td><td>0.27</td><td>0.34</td></tr><tr><td>Mon R1</td><td>11.4 ± 2.2</td><td>0.2 ± 0.2</td><td>9</td><td>1.9</td><td>494</td><td>41</td><td>0.83</td><td>0.17</td><td>0.61</td></tr><tr><td>Rosette</td><td>14.2 ± 3.8</td><td>0.4 ± 0.1</td><td>38</td><td>5.2</td><td>924</td><td>48</td><td>0.75</td><td>0.42</td><td>0.46</td></tr><tr><td>NGC 2264</td><td>15.4 ± 3.2</td><td>0.5 ± 0.2</td><td>44</td><td>5.6</td><td>2091</td><td>42</td><td>0.69</td><td>0.29</td><td>0.38</td></tr><tr><td>M16</td><td>16.6 ± 3.8</td><td>0.5 ± 0.2</td><td>30</td><td>20.1</td><td>6207</td><td>43</td><td>0.49</td><td>0.42</td><td>0.70</td></tr><tr><td>M17</td><td>23.0 ± 6.5</td><td>0.7 ± 0.3</td><td>54</td><td>26.1</td><td>20439</td><td>38</td><td>0.47</td><td>0.42</td><td>0.87</td></tr><tr><td>W48</td><td>14.8 ± 3.1</td><td>0.5 ± 0.3</td><td>91</td><td>137.6</td><td>55013</td><td>100</td><td>0.48</td><td>0.32</td><td>0.89</td></tr><tr><td>Cygnus X South</td><td>14.4 ± 3.6</td><td>0.4 ± 0.2</td><td>75</td><td>19.7</td><td>4744</td><td>119</td><td>0.35</td><td>0.36</td><td>0.66</td></tr><tr><td>Cygnus X North</td><td>17.1 ± 4.3</td><td>0.4 ± 0.2</td><td>130</td><td>34.7</td><td>11383</td><td>194</td><td>0.47</td><td>0.35</td><td>0.64</td></tr><tr><td>NGC 7538</td><td>16.7 ± 4.0</td><td>0.6 ± 0.3</td><td>66</td><td>59.2</td><td>17393</td><td>73</td><td>0.56</td><td>0.40</td><td>0.58</td></tr></tbody></table>

**Note.** Columns show the following: (1) cloud name, (2, 3) median kinetic gas temperature and velocity dispersion for all reliably fit ammonia pixels, (4) number of protostars identified where NH₃ (1,1) integrated intensity is greater than 1 K km s⁻¹, (5) area of map where NH₃ (1,1) integrated intensity is greater than 1 K km s⁻¹, (6) total mass where NH₃ (1,1) integrated intensity is greater than 1 K km s⁻¹, (7) total number of ammonia leaves identified in cloud, (8) fraction of total leaves that are on-filament, (9) fraction of total leaves that are protostellar, (10) fraction of total leaves that are bound ($\alpha_{vir} < 2$).

around $1 \times 10^{22} \text{ cm}^{-2}$, with a minimum of $4.2 \times 10^{21} \text{ cm}^{-2}$ measured in W3-west and a maximum of $2.3 \times 10^{22} \text{ cm}^{-2}$ 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
<figure><img src="image_37.png" /><figcaption>Figure 35. Same as Figure 32 for Mon R1.</figcaption></figure>

<table><caption>Table 7
Leaf Population Statistics</caption><thead><tr><td>Leaf Statistic</td><td>Fraction</td></tr></thead><tbody><tr><td>Bound</td><td>523 of 835 (~63%)</td></tr><tr><td>Starless</td><td>547 of 835 (~66%)</td></tr><tr><td>Protostellar</td><td>288 of 835 (~34%)</td></tr><tr><td>On-filament</td><td>454 of 835 (~54%)</td></tr><tr><td>Bound starless</td><td>294 of 547 (~54%)</td></tr><tr><td>Bound protostellar</td><td>229 of 288 (~80%)</td></tr><tr><td>Bound on-filament</td><td>294 of 454 (~65%)</td></tr><tr><td>Sub-virial (cloud weight pressure)</td><td>573 of 835 (~69%)</td></tr></tbody></table>

Note. Leaf population statistics quoted throughout the paper. All fractions are in relation to the sample of leaves with mass estimates.

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 ~|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.

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.
<figure><img src="image_38.png" /><figcaption>Figure 36. Same as Figure 32 for Rosette.</figcaption></figure>

Errorbars for the total protostars, protostellar surface density, and NH₃ leaves represent the √*N* counting uncertainty. Errorbars for the median kinetic temperature represent the median absolute deviation of each cloud's *T*<sub>K</sub> 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.
<figure><img src="image_39.png" /><figcaption>Figure 37. Same as Figure 32 for NGC 2264.</figcaption></figure>

# 4. Discussion

## 4.1. Leaf/Filament Relationship

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 $\alpha_{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
<figure><img src="image_40.png" /><figcaption>Figure 38. Same as Figure 32 for M16.</figcaption></figure>

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.

We also note the existence of a group of ammonia-identified leaves with uncharacteristically larger masses (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)¹.³³. 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(Mobs) = 3.2 ± 0.5), larger radii (median Reff = 0.7 ± 0.3), warmer temperatures (median TK,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).

## 4.2. Leaf/Protostar Relationship

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
<figure><img src="image_41.png" /><figcaption>**Figure 39.** Same as Figure 32 for M17.</figcaption></figure>

population clearly has larger masses and lower virial parameters than the starless core population. For the 288 protostellar leaves identified, 229 (~80%) have $\alpha_{vir} < 2$. In comparison, 294 of the 547 starless leaves identified (~54%) have $\alpha_{vir} < 2$. As shown in Figure 46, the protostellar population's mass distribution peaks at higher values (median $log(mass) = 1.8 \pm 0.6 M_\odot$) than that of the starless population (median $log(mass) = 1.3 \pm 0.5 M_\odot$), 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 underestimated. This attribute highlights the exceptional environment hubs and ridges provide for cluster formation.

### 4.3. Virial Stability in Low- and High-mass Star-forming Regions

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
<figure><img src="image_42.png" /><figcaption>Figure 40. Same as Figure 32 for W48.</figcaption></figure>

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 ($^{13}$CO, NH$_3$, N$_2$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). $\alpha_{vir}$ ranges from $\sim 10^{-1}$ to $10^2$ and roughly half of the sources fall below the $\alpha_{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 $\alpha_{vir} < 2$. Serpens South and Cepheus-L1251 are closer ($d \sim 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 $\alpha_{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_2H^+$ (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 ($\sim 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 ($\alpha_{vir} < 2$) for clumps. For instance, Svoboda et al. (2016) combined NH$_3$ 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 $\alpha_{vir} < 2$. Similarly, Dunham et al. (2011) found that a separate sample of 456 BGPS clumps had a median $\alpha_{vir} = 0.74$. We note, however, that the NH$_3$ 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$_3$ emission probed by KEYSTONE and included in our leaf catalog. These low-brightness NH$_3$ sources comprise the lower-mass leaves in our sample that have $\alpha_{vir} > 2$ and may be the reason we detect unbound sources that are lacking in the BGPS data.
<figure><img src="image_43.png" /><figcaption>Figure 41. Same as Figure 34 for Cygnus X South.</figcaption></figure>

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 < α<sub>vir</sub> < 1, with a median α<sub>vir</sub> of 0.3. A similar virial analysis of 213 Hi-GAL clumps by Traficante et al. (2018) found that 76% have α<sub>vir</sub> < 1. Using a sample of 1244 ATLASGAL clump detections with masses much larger (typically >10<sup>3</sup>M<sub>☉</sub>) than the leaves presented in this paper, Contreras et al. (2017) found a median α<sub>vir</sub> 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.

### 4.4. Cloud Weight Pressure

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 (Ω<sub>Pw</sub>) 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:

$$ \Omega_{Pw} = -4\pi P_w R^3 \quad (4) $$

$$ \Omega_G = \frac{-1}{2\sqrt{\pi}} \frac{GM^2}{R} \quad (5) $$

$$ \Omega_K = \frac{3}{2} M \sigma^2 \quad (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<sub>w</sub> is cloud weight pressure:

$$ P_w = \pi G \bar{N} N (\mu_H m_H)^2 \quad (7) $$
<figure><img src="image_44.png" /><figcaption>**Figure 42.** Same as Figure 34 for Cygnus X North.</figcaption></figure>

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.<fnref n="31" /> is used to filter spatially each column density map to spatial scales larger than $2^n$ 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.

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 ($\Omega_{Pw}$), kinetic energy ($\Omega_K$), and gravitational potential energy ($\Omega_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

<footnote n="31">³¹ The *atrous.pro* IDL script developed by Erik Rosolowsky, available at https://github.com/low-sky/idl-low-sky/blob/master/wavelet/atrous.pro, was used for this analysis.</footnote>
<figure><img src="image_45.png" /><figcaption>Figure 43. Same as Figure 32 for NGC 7538.</figcaption></figure>

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.

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.

### 4.5. Turbulent Pressure

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 ($\Omega_{\text{Pl}}$) 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
<figure><img src="image_46.png" /><figcaption>**Figure 44.** Heatmap of Pearson correlation coefficients for the 10 variables of interest examined: leaf on-filament fraction, “bound” leaf fraction (i.e., the fraction of leaves with $\alpha_{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 colorbar ranges from -1 (perfect anti-correlation) to 1 (perfect positive correlation). Panels with a black outline denote statistically significant correlations that are displayed in Figure 45.</figcaption></figure>

pressure). $\Omega_{Pt}$ is calculated from Equation (4), with $P_w$ being replaced with $P_t$ given by

$$P_T = \mu_H m_H \times \rho_{C^{18}O} \times \sigma_{C^{18}O}^2, \quad (8)$$

where $\sigma_{C^{18}O}$ is the velocity dispersion measured from C$^{18}$O (3–2), a moderate density tracer, and $\rho_{C^{18}O}$ is the volume density at which the C$^{18}$O (3–2) emission originates. Here, we assume the C$^{18}$O (3–2) emission is tracing a volume density of $3 \times 10^4 \text{ cm}^{-3}$, which is the critical density ($n_{cr}(u-l) = A_{ul}/\gamma_{ul}$, where $A_{ul}$ is the Einstein A coefficient and $\gamma_{ul}$ is the collisional rate coefficient for collisions with H$_2$) of C$^{18}$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$^{18}$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$^{18}$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$^{-1}$ for the velocity dispersion. After the line fitting, pixels must meet the following criteria to be included in our final parameter maps:

1. $\sigma > 0.05 \text{ km s}^{-1}$ (below $0.05 \text{ km s}^{-1}$ is unrealistic since our channel width is only $\sim 0.11 \text{ km s}^{-1}$;

2. $T_{\text{peak, err}} < 1 \text{ K}$;

3. $\sigma_{\text{err}} < 0.5 \text{ km s}^{-1}$;

4. $V_{\text{LSR, err}} < 0.75 \text{ km s}^{-1}$.

The final parameter maps for the velocity dispersion, $\sigma_{C^{18}O}$, are shown in Figures 49 and 50 and . The 52 leaves with at least one reliably fit C$^{18}$O (3–2) pixel are shown by red contours in Figures 49 and 50 and . The median value of $\sigma_{C^{18}O}$ is calculated within the boundaries of each leaf and converted into a turbulent pressure using Equations (4) and (8). Although most of the C$^{18}$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
<figure><img src="image_47.png" /></figure>

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

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

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

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

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

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

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

**Figure 45.** Statistically significant correlations found from Figure 44 when neglecting data errorbars. The dense gas mass, number of protostars, and surface protostar density are calculated within the boundaries of the KEYSTONE-mapped regions in each cloud where NH₃ (1,1) integrated intensity was greater than 1.0 K km s⁻¹. The errorbars on each data point are calculated as described in Section 3.7.
<figure><img src="image_55.png" /><figcaption>**Figure 46.** Stacked histograms of mass (upper left), effective radius (upper right), average kinetic temperature (lower left), and average velocity dispersion (lower right) for the leaves identified as on-filament, off-filament, protostellar, starless, and hub/ridge. The median and median absolute deviation of the distributions are shown in the upper right corner of each panel. The black vertical lines show the distribution median.</figcaption></figure>

as sharp increases in $\sigma_{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.

The right panel of Figure 48 shows the virial plane for structures with C¹⁸O (3–2) measurements when using $\Omega_{Pl}$ as the pressure term in the virial equation. All of these structures fall within the pressure-dominated, sub-virial quadrant. The ratio of $\Omega_{Pl}/\Omega_{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_{C18O}$, we also calculate $P_T$ using a factor of 10 lower value ($3 \times 10^3 \text{ cm}^{-3}$) for $\rho_{C18O}$ 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 $\Omega_{Pl}/\Omega_{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_{C18O}$ 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
<figure><img src="image_56.png" /><figcaption>**Figure 47.** Left: H₂ column density map of NGC 7538 with white contours at 4 × 10²¹ cm⁻² and 8 × 10²¹ cm⁻². Black outlines the extent of the KEYSTONE mapping of the region. Right: spatially filtered H₂ column density map over the cyan dotted area overlaid in the left panel. The map includes spatial scales larger than 16 pixels, which is equivalent to ~96'' or ~1.3 pc at the distance of NGC 7538.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>**Figure 48.** Left: virial plane for all ammonia-identified leaves displayed in Figure 31 showing the balance between three energy densities in the virial equation: gravitational (Ω<sub>G</sub), cloud weight pressure (Ω<sub>Pw</sub>), and kinetic (Ω<sub>K</sub>). Sources to the right of the vertical line are sub-virial, while sources to the left are super-virial. Sources above the horizontal line are gravitationally dominated, while those below the line are pressure-dominated. Right: virial plane using instead the turbulent pressure energy density (Ω<sub>Pt</sub>) for all ammonia-identified leaves with reliable C¹⁸O (3–2) velocity dispersion measurements. The data points are calculated using the critical density (3 × 10⁷ cm⁻³) of C¹⁸O (3–2) as the volume density traced by its emission. The errorbars show where the source would fall on the plot if the density traced by the C¹⁸O were a factor of 10 lower (3 × 10³ cm⁻³).</figcaption></figure>

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.

**5. Summary**

We have presented initial observations and results from KFPA Examinations of Young STellar Object Natal Environments (KEYSTONE), a survey of filamentary dense gas structure in 11 GMCs (Cygnus X North, Cygnus X South, M16, M17, Mon R1, Mon R2, NGC 2264, NGC 7538, Rosette, W3, and W48) at distances of 0.9–3.0 kpc. We identified 856 dense gas clumps, traced by NH₃ (1,1) emission, across all the observed clouds using a dendrogram analysis. Simultaneous line fitting of the NH₃ (1,1) and (2,2) emission provided estimates of kinetic gas temperature, centroid velocity, velocity dispersion, and para-NH₃ column density for the dense gas. These parameter maps and the NH₃ (1,1) and (2,2) data cubes are publicly available (see footnote 26).

The ammonia parameter maps were used to derive virial stability parameters for each dense gas structure identified, providing insight into whether or not the gravitational potential energy of the structures is enough to overcome their internal
<figure><img src="image_58.png" /><figcaption>**Figure 49.** C<sup>18</sup>O (3–2) velocity dispersion measured from Gaussian fits to JCMT observations of DR21 in Cygnus X North (top left), G79.34 in Cygnus X South (top right), W3(OH) (bottom left), and W3-Main (bottom right). Red contours denote ammonia-identified leaves that have at least one reliably fit C<sup>18</sup>O (3–2) pixel and were included in our turbulent pressure virial analysis. Blue contours show all other ammonia-identified leaves in the field of view.</figcaption></figure>

kinetic energies in the absence of magnetic fields or external pressure. HOBYS Herschel observations of dust continuum emission were utilized to create H<sub>2</sub> column density maps, identify young protostellar candidates, and identify filamentary structures in each region. JCMT observations of C<sup>18</sup>O (3–2) emission were accessed to determine the turbulent pressure applied by the ambient molecular cloud upon a subset of the dendrogram-identified dense gas structures. Our main results are listed below.

1. Significant variations in kinetic gas temperature are observed between clouds, with median T<sub>K</sub> = 11.4 ± 2.2 K in the coldest region (Mon R1) and T<sub>K</sub> = 23.0 ± 6.5 K in the warmest (M17). The velocity dispersion distributions are more similar between clouds, with characteristic median values of 0.3–0.7 km s<sup>-1</sup>.

2. Of the 835 ammonia-identified clumps with mass estimates, 523 (~63%) have virial parameters less than

two, suggesting the mass of those structures is gravitationally bound and more susceptible to gravitational collapse when neglecting the effects of magnetic fields or external pressure. Similar analyses in nearby low-mass star-forming clouds have found much higher rates of gravitationally bound ammonia-identified cores, which may suggest ammonia is more widespread in GMCs than in nearby clouds or that gravity is more important to structure stability at small scales.

3. The fraction of ammonia-identified clumps that are spatially coincident with filaments identified in the H<sub>2</sub> column density maps ranges from 0.35 in Cygnus X South to 1.0 in W3-west. These values are consistent with core on-filament fractions found from dust continuum observations of nearby star-forming regions, which tend to be from ~0.4 to 0.8 depending on the cloud and core class considered.

4. On- and off-filament clumps show no substantial differences in their virial parameter, mass, radius,
<figure><img src="image_59.png" /><figcaption>**Figure 50.** Same as Figure 49 for M16 (top left), M17 (top right), and NGC 7538 (bottom).</figcaption></figure>

temperature, and velocity dispersion distributions. We do find, however, a tendency for clouds with low dense gas mass to have a higher fraction of on-filament clumps. These findings may indicate that filaments play a lesser role in the star formation process of high-mass GMCs. In those environments, dense gas may be more widespread, allowing for clump formation to be equally as likely on and off filaments. In lower-mass environments where dense gas is less widespread, however, clump formation may be limited to the filaments that harbor the main supply of dense gas.

5. In several regions there are “hubs” or “ridges” of dense gas that have much higher masses and lower virial parameters than the other clumps in their respective cloud. These hubs and ridges tend to be located at the intersections of multiple filaments or located near/within a single filament, are often associated with H₂O maser emission, and typically host multiple protostars. Based on these characteristics, hubs may be the sites of future cluster formation.

6. When considering the external pressure exerted on the clumps, most are considered sub-virial and pressure-dominated structures. This characteristic state may indicate that high-mass clumps spend the majority of their lifetime confined by external pressure. Over time, as the clumps accrete mass from their surroundings, they may gain enough mass to be gravitationally dominated and undergo gravitational collapse or fragmentation.

## 6. Future Work

Although it was not the focus of this paper, a key use-case of the KEYSTONE data is the analysis of filament kinematics in GMCs. For instance, with regard to the observed spatial relationship between MYSOs and stellar clusters with filament intersections, the KEYSTONE data could be used to determine: (1) whether or not the observed clumps and filaments are truly velocity-coherent structures (Pineda et al. 2010; Chen et al. 2019a), and (2) whether the mass flow rates along them are large enough to produce MYSOs. Several independent studies
have already measured gas motions in individual filaments such as Serpens South (Friesen et al. 2013; Kirk et al. 2013a), G035.39-00.33 (Henshaw et al. 2013), DR21 (Schneider et al. 2010a), W43-MM1 (Nguyen-Lu'o'ng et al. 2013), M17SW (Chen et al. 2019b), and eight other high-mass filaments (Lu et al. 2018), noting that the observed mass flow rates could supply the mass required to assemble the stellar clusters at their centers. Similarly, observations of the Large Magellanic Cloud by Fukui et al. (2015) showed two separate instances of MYSOs forming at the centers of adjoining filaments that have gas flowing into the central junction. Svoboda et al. (2016) and Motte et al. (2018a) also contend that such mass flow onto sites of cluster formation is prominent throughout Galactic high-mass star-forming regions, providing the mass build-up and compression necessary to form stellar clusters. These studies, however, have focused primarily on regions that have already formed clusters/MYSOs and do not address the conditions of the parental clumps, i.e., the dense gas out of which stellar clusters form, prior to the onset of star formation. As such, the gas velocity patterns of those regions are susceptible to distortions due to stellar feedback, which raises questions about the applicability of their kinematic measurements. Furthermore, observations after the onset of star formation cannot be used to decipher whether clumps form before or after filaments collide to form intersections. Using the clump catalog presented in this paper, in conjunction with an analysis of the adjacent filament kinematics, it is now possible to investigate these questions across multiple high-mass star-forming environments.

In addition to dense gas kinematics, the KEYSTONE data can provide insight into temperature and chemistry variations as a function of environment within GMCs. For example, the KEYSTONE observations of the NH₃ (1,1), (2,2), (3,3), (4,4), and (5,5) transitions probe a range of gas temperatures up to ~200 K. Since dust temperatures derived from Herschel data become increasingly uncertain above 20 K as the SED peak moves toward wavelengths shorter than 160 µm (Chen et al. 2016), ammonia-derived gas temperatures will be important for understanding how the OB associations are impacting the star formation in the KEYSTONE target clouds. Furthermore, the KEYSTONE observations will constrain the abundances of NH₃ in GMCs. These abundances, in combination with the gas temperatures, will provide a way to understand how temperature impacts the formation/destruction of ammonia in GMCs.

J.K., J.D.F., E.R., and M.C.C. acknowledge the financial support of a Discovery Grant from NSERC of Canada. S.B. and N.S. acknowledge support by the French ANR and the German DFG through the project "GENESIS" (ANR-16-CE92-0035-01/DFG1591/2-1). The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating

universities in the United Kingdom and Canada. This research made use of astrodendro, a Python package to compute dendrograms of Astronomical data (http://www.dendrograms.org/), Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), and pyspeckit (https://pyspeckit.readthedocs.io/en/latest/), a Python spectroscopic analysis and reduction toolkit (Ginsburg & Mirocha 2011).

Facilities: GBT, Herschel, JCMT.

## Appendix

### Distance Dependence of Virial Parameters

The virial analysis presented in this paper used the native KEYSTONE resolution for all clouds analyzed, despite the cloud distances ranging from 0.9 to 3 kpc. Such distance variations provide a factor of ~3 range of linear spatial resolutions, which may lead to different types of structures (in terms of size and mass) being identified in each cloud. To test whether or not the distance dependence has any influence on the main results of this paper, we convolved the NH₃ (1,1) and (2,2) cubes for NGC 2264, Mon R1, and Mon R2 ($d = 0.9$ kpc, $\theta \sim 0.13$ pc), the closest clouds in KEYSTONE, to the linear resolution of W48 ($d = 3.0$ kpc, $\theta \sim 0.45$ pc), the most distant cloud observed. This process degrades the resolution of the observations by a factor of 3.0/0.9 ~ 3.3 to a final resolution of ~103". Following a similar distance bias analysis by Baldeschi et al. (2017) on Herschel maps, we also down-sampled the cubes by the same factor along each spatial axis. Finally, Gaussian noise with a mean of zero and standard deviation of $s_N\sqrt{1-d_0/d_1}$ was added to the cubes, where $s_N$ is the median of the original cube's rms map, $d_0$ is 0.9 pc and $d_1$ is 3.0 kpc. This noise addition attempts to return the native rms level to the convolved and resampled cubes, which have lower noise levels due to the spatial averaging.

These processes allowed us to simulate what might be observed if NGC 2264, Mon R1, and Mon R2 were at a distance of 3 kpc. The convolved, downsampled, noise-injected cubes were then run through the line-fitting pipeline to produce a new set of ammonia parameter maps and integrated intensity maps. Finally, we repeated the virial analysis presented in Section 3.5 using the new set of maps and assuming their distance is 3.0 kpc.

In Figures 51–53, we compare the results of the original virial analyses in NGC 2264, Mon R1, and Mon R2 to their distance-adjusted virial analysis. Leaves from the original analysis that fall within a leaf identified in the distance-adjusted analysis are tagged with a specific color in each plot of Figures 51–53. As expected, much fewer structures are identified by the dendrogram in the distance-adjusted analysis. Many of the small structures in the original analysis are lumped together into a single large structure in the distance-adjusted analysis due to the lower resolution. In terms of virial parameters, the structures in the distance-adjusted analysis are all gravitationally bound. This distance bias is likely why more bound structures are observed in W48 than these closer clouds. At 900 pc, NGC 2264, Mon R1, and Mon R2 represent the extreme examples of the distance bias in our sample. The impact will undoubtedly be less for the moderate-distance clouds in our sample (e.g., Cygnus X, W3, NGC 7538, etc.), for which the distance difference is lower. This statement is echoed by the fact that many of the cloud attributes analyzed in this paper (e.g., dense gas mass, leaf on-filament fraction,
<figure><img src="image_60.png" /><figcaption>**Figure 51.** Top row: leaf virial parameters in NGC 2264 using native KEYSTONE resolution (left) and after convolving the NH₃ (1,1) and (2,2) cubes to the linear resolution of W48 (θ ~ 0.45 pc) (right), downsampling the number of pixels, adding white noise, and re-running the full analysis. Bottom: integrated intensity map obtained using the convolved NH₃ (1,1) cube. Ellipses represent the peaks of leaves identified when using the native resolution data. The green dendrogram masks show the extent of the leaves identified using the convolved data. Ellipses falling within a dendrogram mask are tagged as a colored pair in the top row plots.</figcaption></figure>
<figure><img src="image_61.png" /><figcaption>Figure 52. Same as Figure 51 for Mon R1.</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>Figure 53. Same as Figure 51 for Mon R2.</figcaption></figure>

number of leaves, etc.) do not depend strongly on cloud distance (see Section 3.7 and Figure 44). Nevertheless, we are indeed tracing different scales in the star formation hierarchy by including a range of cloud distances in our analysis. Future high-resolution observations of the more distant KEYSTONE targets with the Karl G. Jansky Very Large Array (VLA) or Next Generation VLA (ngVLA) could provide the means to compare cloud structures between clouds on more similar spatial scales.

## ORCID iDs

Jared Keown https://orcid.org/0000-0003-2628-0250

James Di Francesco https://orcid.org/0000-0002-9289-2450

Erik Rosolowsky https://orcid.org/0000-0002-5204-2259

Ayushi Singh https://orcid.org/0000-0001-6022-3618

Helen Kirk https://orcid.org/0000-0002-5779-8549

L. D. Anderson https://orcid.org/0000-0001-8800-1793

Michael Chun-Yuan Chen https://orcid.org/0000-0003-4242-973X

Davide Elia https://orcid.org/0000-0002-9120-5890

Rachel Friesen https://orcid.org/0000-0001-7594-8128

Adam Ginsburg https://orcid.org/0000-0001-6431-9633

Stefano Pezzuto https://orcid.org/0000-0001-7852-1971

Paola Caselli https://orcid.org/0000-0003-1481-7911

Hong-Li Liu https://orcid.org/0000-0003-3343-9645

Steven Longmore https://orcid.org/0000-0001-6353-0170

Philip C. Myers https://orcid.org/0000-0002-2885-1806

Stella S. R. Offner https://orcid.org/0000-0003-1252-9916

Patricio Sanhueza https://orcid.org/0000-0002-7125-7685

Nicola Schneider https://orcid.org/0000-0003-3485-6678

Ian Stephens https://orcid.org/0000-0003-3017-4418

## References

Adams, M. T., Strom, K. M., & Strom, S. E. 1979, ApJ, 230, L183

André, P., Men'schikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Auddy, S., Myers, P. C., Basu, S., et al. 2019, ApJ, 872, 207

Baldeschi, A., Elia, D., Molinari, S., et al. 2017, MNRAS, 466, 3682

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

Baxter, E. J., Covey, K. R., Muench, A. A., et al. 2009, AJ, 138, 963

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

Billington, S. J., Urquhart, J. S., Figura, C., Eden, D. J., & Moore, T. J. T. 2019, MNRAS, 483, 3146

Blitz, L., & Stark, A. A. 1986, ApJ, 300, L89

Blitz, L., & Thaddeus, P. 1980, ApJ, 241, 676

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bonatto, C., Santos, J. F. C., Jr., & Bica, E. 2006, A&A, 445, 567

Bonnell, I. A., Bate, M. R., & Vine, S. G. 2003, MNRAS, 343, 413

Bontemps, S., Motte, F., Csengeri, T., & Schneider, N. 2010, A&A, 524, A18

Boyden, R. D., Koch, E. W., Rosolowsky, E. W., & Offner, S. S. R. 2016, ApJ, 833, 233
Boyden, R. D., Offner, S. S. R., Koch, E. W., & Rosolowsky, E. W. 2018, ApJ, 860, 157

Bresnahan, D., Ward-Thompson, D., Kirk, J. M., et al. 2018, A&A, 615, A125

Burkhart, B. 2018, ApJ, 863, 118

Campbell, B. 1984, ApJL, 282, L27

Carpenter, J. M., Meyer, M. R., Dougados, C., Strom, S. E., & Hillenbrand, L. A. 1997, AJ, 114, 198

Castelaz, M. W., Grasdalen, G. L., Hackwell, J. A., Capps, R. W., & Thompson, D. 1985, AJ, 90, 1113

Chen, H. H.-H., Pineda, J. E., Goodman, A. A., et al. 2019a, ApJ, 877, 93

Chen, H.-R. V., Zhang, Q., Wright, M. C. H., et al. 2019b, ApJ, 875, 24

Chen, M. C.-Y., Di Francesco, J., Johnstone, D., et al. 2016, ApJ, 826, 95

Ching, T.-C., Lai, S.-P., Zhang, Q., et al. 2017, ApJ, 838, 121

Colley, D. 1980, MNRAS, 193, 495

Comerón, F., Djupvik, A. A., Torra, J., Schneider, N., & Pasquali, A. 2019, A&A, 622, A134

Comerón, F., & Pasquali, A. 2012, A&A, 543, A101

Contreras, Y., Rathborne, J. M., Guzman, A., et al. 2017, MNRAS, 466, 340

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

Csengeri, T., Bontemps, S., Schneider, N., Motte, F., & Dib, S. 2011b, A&A, 527, A135

Cunningham, A. J., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2018, MNRAS, 476, 771

Dahm, S. E. 2008, in Handbook of Star Forming Regions Vol. 1, ed. B. Reipurth (San Francisco, CA: ASP), 966

Didelon, P., Motte, F., Tremblin, P., et al. 2015, A&A, 584, A4

Di Francesco, J., Evans, N. J., II, Caselli, P., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 17

Di Francesco, J., Sadavoy, S., Motte, F., et al. 2010, A&A, 518, L91

Dobashi, K., Matsumoto, T., Shimoikura, T., et al. 2014, ApJ, 797, 58

Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2013, A&A, 558, A125

Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2014, A&A, 570, A1

Dunham, M. K., Rosolowsky, E., Evans, N. J., II, et al. 2010, ApJ, 717, 1157

Dunham, M. K., Rosolowsky, E., Evans, N. J., II, Cyganowski, C., & Urquhart, J. S. 2011, ApJ, 741, 110

Elia, D., & Pezzuto, S. 2016, MNRAS, 461, 1328

Elmegreen, B. G., Lada, C. J., & Dickinson, D. F. 1979, ApJ, 230, 415

Evans, C. J., Smartt, S. J., Lee, J.-K., et al. 2005, A&A, 437, 467

Falgarone, E., & Puget, J. L. 1986, A&A, 162, 235

Fallscheer, C., Reid, M. A., Di Francesco, J., et al. 2013, ApJ, 773, 102

Federrath, C. 2016, MNRAS, 457, 375

Field, G. B., Blackman, E. G., & Keto, E. R. 2011, MNRAS, 416, 710

Friesen, R. K., Bourke, T. L., Di Francesco, J., Guttermuth, R., & Myers, P. C. 2016, ApJ, 833, 204

Friesen, R. K., Medeiros, L., Schnee, S., et al. 2013, MNRAS, 436, 1513

Friesen, R. K., Pineda, J. E., Rosolowsky, E., et al. 2017, ApJ, 843, 63

Fukui, Y., Harada, R., Tokuda, K., et al. 2015, ApJL, 807, L4

Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

Gatley, I., Becklin, E. E., Sellgren, K., & Werner, M. W. 1979, ApJ, 233, 575

Ginsburg, A., Glenn, J., Rosolowsky, E., et al. 2013, ApJS, 208, 14

Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit v0.1.19, Astrophysics Source Code Library, ascl:1109.001

Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., et al. 2009, Natur, 457, 63

Gritschneder, M., Burkert, A., Naab, T., & Walch, S. 2010, ApJ, 723, 971

Guesten, R., & Fiebig, D. 1988, A&A, 204, 253

Guilloteau, S., Wilson, T. L., Martin, R. N., Batria, W., & Pauls, T. A. 1983, A&A, 124, 322

Hachisuka, K., Brunthaler, A., Menten, K. M., et al. 2006, ApJ, 645, 337

Hansen, C. E., Klein, R. I., McKee, C. F., & Fisher, R. T. 2012, ApJ, 747, 22

Hanson, M. M. 2003, ApJ, 597, 957

Hennebelle, P. 2013, A&A, 556, A153

Hennemann, M., Motte, F., Schneider, N., et al. 2012, A&A, 543, L3

Henshaw, J. D., Caselli, P., Fontani, F., et al. 2013, MNRAS, 428, 3425

Hester, J. J., Scowen, P. A., Sankrit, R., et al. 1996, AJ, 111, 2349

Heyer, M. H., Williams, J. P., & Brunt, C. M. 2006, ApJ, 643, 956

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

Hill, T., André, P., Arzoumanian, D., et al. 2012a, A&A, 548, L6

Hill, T., Motte, F., Didelon, P., et al. 2012b, A&A, 542, A114

Ho, P. T. P., Barrett, A. H., Myers, P. C., et al. 1979, ApJ, 234, 912

Ho, P. T. P., & Townes, C. H. 1983, ARA&A, 21, 239

Hoffman, I. M. 2012, ApJ, 759, 76

Hoffman, I. M., & Seojin Kim, S. 2011, ApJL, 739, L15

Hoffmeister, V. H., Chini, R., Scheyda, C. M., et al. 2008, ApJ, 686, 310

Hogge, T., Jackson, J., Stephens, I., et al. 2018, ApJS, 237, 27

Johnston, K. J., Stolovy, S. R., Wilson, T. L., Henkel, C., & Mauersberger, R. 1989, ApJL, 343, L41

Kallas, E., & Reich, W. 1980, A&AS, 42, 227

Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., II, & Lee, C. W. 2008, A&A, 487, 993

Kauffmann, J., & Pillai, T. 2010, ApJL, 723, L7

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

Keown, J., Di Francesco, J., Kirk, H., et al. 2017, ApJ, 850, 3

Kerr, R., Kirk, H., Di Francesco, J., et al. 2019, ApJ, 874, 147

Kirk, H., Friesen, R. K., Pineda, J. E., et al. 2017, ApJ, 846, 144

Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009

Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013a, ApJ, 766, 115

Kirk, J. M., Ward-Thompson, D., Palmeirim, P., et al. 2013b, MNRAS, 432, 1424

Knödltseder, J. 2000, A&A, 360, 539

Koch, E. W., Ward, C. G., Offner, S., Loeppky, J. L., & Rosolowsky, E. W. 2017, MNRAS, 471, 1506

Könyves, V., André, P., Men'shchikov, A., et al. 2010, A&A, 518, L106

Könyves, V., André, P., Men'shchikov, A., et al. 2015, A&A, 584, A91

Krumholz, M. R., & McKee, C. F. 2008, Natur, 451, 1082

Kryukova, E., Megeath, S. T., Hora, J. L., et al. 2014, AJ, 148, 11

Kun, M., Kiss, Z. T., & Balog, Z. 2008, in Star-forming Regions in Cepheus, ed. B. Reipurth (San Francisco, CA: ASP), 136

Kurono, Y., Saito, M., Kamazaki, T., Morita, K.-I., & Kawabe, R. 2013, ApJ, 765, 85

Kutner, M. L., Dickman, R. L., Tucker, K. D., & Machnik, D. E. 1979, ApJ, 232, 724

Lada, C. J. 1976, ApJS, 32, 603

Lada, C. J., Forbrich, J., Lombardi, M., & Alves, J. F. 2012, ApJ, 745, 190

Lada, E. A. 1992, ApJL, 393, L25

Lang, K. R., & Willson, R. F. 1980, ApJ, 238, 867

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

Lee, K. I., Fernández-López, M., Storm, S., et al. 2014, ApJ, 797, 76

Li, D., Kauffmann, J., Zhang, Q., & Chen, W. 2013, ApJL, 768, L5

Li, G.-X., Urquhart, J. S., Leurini, S., et al. 2016, A&A, 591, A5

Liu, T., Li, P. S., Juvela, M., et al. 2018, ApJ, 859, 151

Lombardi, M., Alves, J., & Lada, C. J. 2011, A&A, 535, A16

Lovas, F. J. 2004, JPCRD, 33, 177

Lu, X., Zhang, Q., Liu, H. B., et al. 2018, ApJ, 855, 9

Mairs, S., Johnstone, D., Offner, S. S. R., & Schnee, S. 2014, ApJ, 783, 60

Mallick, K. K., Ojha, D. K., Tamura, M., et al. 2014, MNRAS, 443, 3218

Mangum, J. G., Emerson, D. T., & Greisen, E. W. 2007, A&A, 474, 679

Mangum, J. G., Wootten, A., & Mundy, L. G. 1991, ApJ, 378, 576

Margulis, M., & Lada, C. J. 1986, ApJL, 309, L87

Marsh, K. A., Kirk, J. M., André, P., et al. 2016, MNRAS, 459, 342

Mckee, C. F. 1989, ApJ, 345, 782

Megeath, S. T., Herter, T., Beichman, C., et al. 1996, A&A, 307, 775

Megeath, S. T., Townsley, L. K., Oey, M. S., & Tieftrunk, A. R. 2008, in Handbook of Star Forming Regions Vol. 1, ed. B. Reipurth (San Francisco, CA: ASP), 264

Men'shchikov, A. 2013, A&A, 560, A63

Men'shchikov, A., André, P., Didelon, P., et al. 2010, A&A, 518, L103

Men'shchikov, A., André, P., Didelon, P., et al. 2012, A&A, 542, A81

Merello, M., Molinari, S., Rygl, K. L. J., et al. 2019, MNRAS, 483, 5355

Minier, V., Booth, R. S., & Conway, J. E. 2000, A&A, 362, 1093

Minier, V., Tremblin, P., Hill, T., et al. 2013, A&A, 550, A50

Miville-Deschênes, M.-A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104

Molinari, S., Schisano, E., Faustini, F., et al. 2011, A&A, 530, A133

Molinari, S., Swinyard, B., Bally, J., et al. 2010, PASP, 122, 314

Montalban, J., Bachiller, R., Martin-Pintado, J., Tafalla, M., & Gomez-Gonzalez, J. 1990, A&A, 233, 527

Moscadelli, L., Reid, M. J., Menten, K. M., et al. 2009, ApJ, 693, 406

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

Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243

Motte, F., Nony, T., Louvet, F., et al. 2018b, NatAs, 2, 478

Motte, F., Zavagno, A., Bontemps, S., et al. 2010, A&A, 518, L77

Myers, P. C. 2009, ApJ, 700, 1609

Nakano, M., Soejima, T., Chibueze, J. O., et al. 2017, PASJ, 69, 16

Naranjo-Romero, R., Vázquez-Semadeni, E., & Loughnane, R. M. 2015, ApJ, 814, 48

Nguyen Luong, Q., Motte, F., Hennemann, M., et al. 2011, A&A, 535, A76

Nguyen-Lu'o'ng, Q., Motte, F., Carlhoff, P., et al. 2013, ApJ, 775, 88

Nishimura, A., Costes, J., Inaba, T., et al. 2017, PASJ, submitted (arXiv:1706.06002)
Nishimura, A., Minamidani, T., Umemoto, T., et al. 2018, PASJ, 70, S42

Oey, M. S., Watson, A. M., Kern, K., & Walth, G. L. 2005, AJ, 129, 393

Offner, S. S. R., Bisbas, T. G., Viti, S., & Bell, T. A. 2013, ApJ, 770, 49

Offner, S. S. R., & Chaban, J. 2017, ApJ, 847, 104

Offner, S. S. R., & Liu, Y. 2018, NatAs, 2, 896

Ogura, K. 1984, PASJ, 36, 139

Ohashi, S., Sanhueza, P., Chen, H.-R. V., et al. 2016, ApJ, 833, 209

Ojha, D. K., Tamura, M., Nakajima, Y., et al. 2004, ApJ, 608, 797

Oliveira, J. M. 2008, in Handbook of Star Forming Regions Vol. 2, ed. B. Reipurth (San Francisco, CA: ASP), 599

Oliver, R. J., Masheder, M. R. W., & Thaddeus, P. 1996, A&A, 315, 578

Olmi, L., Araya, E. D., Chapin, E. L., et al. 2011, ApJ, 715, 1132

Pagani, L. P., & Nguyen-Q-Rieu 1987, A&A, 181, 112

Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38

Pattle, K., Ward-Thompson, D., Kirk, J. M., et al. 2015, MNRAS, 450, 1094

Pattle, K., Ward-Thompson, D., Kirk, J. M., et al. 2017, MNRAS, 464, 4255

Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979

Pérez-Beaupuits, J. P., Stutzki, J., Ossenkopf, V., et al. 2015, A&A, 575, A9

Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72

Pillai, T., Caselli, P., Kauffmann, J., et al. 2012, ApJ, 751, 135

Pillai, T., Kauffmann, J., Tan, J. C., et al. 2015, ApJ, 799, 74

Pillai, T., Kauffmann, J., Wyrowski, F., et al. 2011, A&A, 530, A118

Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, ApJL, 712, L116

Pirogov, L. E. 2009, ARep, 53, 1127

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11

Pokhrel, R., Gutermuth, R., Ali, B., et al. 2016, MNRAS, 461, 22

Polychroni, D., Schisano, E., Elia, D., et al. 2013, ApJL, 777, L33

Povich, M. S., Churchwell, E., Bieging, J. H., et al. 2009, ApJ, 696, 1278

Purcell, C. R., Longmore, S. N., Walsh, A. J., et al. 2012, MNRAS, 426, 1972

Qiu, K., Zhang, Q., & Menten, K. M. 2011, ApJ, 728, 6

Rapson, V. A., Pipher, J. L., Gutermuth, R. A., et al. 2014, ApJ, 794, 124

Rayner, T. S. M., Griffin, M. J., Schneider, N., et al. 2017, A&A, 607, A22

Redaelli, E., Alves, F. O., Caselli, P., et al. 2017, ApJ, 850, 202

Reid, M. A., & Wilson, C. D. 2005, ApJ, 625, 891

Rivera-Gálvez, S., Román-Zúñiga, C. G., Jiménez-Bailón, E., et al. 2015, AJ, 150, 191

Rivera-Ingraham, A., Martin, P. G., Polychroni, D., et al. 2013, ApJ, 766, 85

Rivera-Ingraham, A., Martin, P. G., Polychroni, D., et al. 2015, ApJ, 809, 81

Rivera-Ingraham, A., Martin, P. G., Polychroni, D., & Moore, T. J. T. 2011, ApJ, 743, 39

Rivera-Ingraham, A., Ristorcelli, I., Juvela, M., et al. 2016, A&A, 591, A90

Román-Zúñiga, C. G., & Lada, E. A. 2008, in Handbook of Star Forming Regions Vol. 1, ed. B. Reipurth (San Francisco, CA: ASP), 928

Rosolowsky, E., Dunham, M. K., Ginsburg, A., et al. 2010, ApJS, 188, 123

Rosolowsky, E., & Leroy, A. 2006, PASP, 118, 590

Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338

Russeil, D., Figureira, M., Zavagno, A., et al. 2019, A&A, 625, 134

Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79

Rygl, K. L. J., Goedhart, S., Polychroni, D., et al. 2014, MNRAS, 440, 427

Rygl, K. L. J., Wyrowski, F., Schuller, F., & Menten, K. M. 2010, A&A, 515, A42

Sandell, G., Goss, W. M., & Wright, M. 2005, ApJ, 621, 839

Sandell, G., Wright, M., & Forster, J. R. 2003, ApJL, 590, L45

Schilke, P., Mauersberger, R., Walmsley, C. M., & Wilson, T. L. 1990, A&A, 227, 220

Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29

Schneider, N., Bontemps, S., Girichidis, P., et al. 2015, MNRAS, 453, L41

Schneider, N., Bontemps, S., Motte, F., et al. 2016, A&A, 591, A40

Schneider, N., Bontemps, S., Simon, R., et al. 2006, A&A, 458, 855

Schneider, N., Bontemps, S., Simon, R., et al. 2011, A&A, 529, A1

Schneider, N., Csengeri, T., Bontemps, S., et al. 2010a, A&A, 520, A49

Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11

Schneider, N., Motte, F., Bontemps, S., et al. 2010b, A&A, 518, L83

Schneider, N., Stutzki, J., Winnewisser, G., & Block, D. 1998, A&A, 335, 1049

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

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

Scoville, N. Z., Sargent, A. I., Sanders, D. B., et al. 1986, ApJ, 303, 416

Seifried, D., & Walch, S. 2015, MNRAS, 452, 2410

Seo, Y. M., Shirley, Y. L., Goldsmith, P., et al. 2015, ApJ, 805, 185

Shirley, Y. L. 2015, PASP, 127, 299

Slysh, V. I., Kalenskii, S. V., & Val'tts, I. E. 1995, ApJ, 442, 668

Smith, R. J., Glover, S. C. O., Clark, P. C., Klessen, R. S., & Springel, V. 2014a, MNRAS, 441, 1628

Smith, R. J., Glover, S. C. O., & Klessen, R. S. 2014b, MNRAS, 445, 2900

Stephens, I. W., Jackson, J. M., Whitaker, J. S., et al. 2016, ApJ, 824, 29

Stutzki, J., & Guesten, R. 1990, ApJ, 356, 513

Stutzki, J., Stacey, G. J., Genzel, R., et al. 1988, ApJ, 332, 379

Sugiyama, K., Fujisawa, K., Doi, A., et al. 2008, PASJ, 60, 23

Surcis, G., Vlemmings, W. H. T., van Langevelde, H. J., & Hutawarakorn Kramer, B. 2012, A&A, 541, A47

Svoboda, B. E., Shirley, Y. L., Battersby, C., et al. 2016, ApJ, 822, 59

Tan, J. C., Kong, S., Butler, M. J., Caselli, P., & Fontani, F. 2013, ApJ, 779, 96

Teixeira, P. S., Lada, C. J., Marengo, M., & Lada, E. A. 2012, A&A, 540, A83

Thronson, H. A., Jr., & Lada, C. J. 1983, ApJ, 269, 175

Tieftrunk, A. R., Gaume, R. A., Claussen, M. J., Wilson, T. L., & Johnston, K. J. 1997, A&A, 318, 931

Tieftrunk, A. R., Megeath, S. T., Wilson, T. L., & Rayner, J. T. 1998, A&A, 336, 991

Tigé, J., Motte, F., Russeil, D., et al. 2017, A&A, 602, A77

Traficante, A., Duarte-Cabral, A., Elia, D., et al. 2018, MNRAS, 477, 2220

Tremblin, P., Minier, V., Schneider, N., et al. 2013, A&A, 560, A19

Tremblin, P., Schneider, N., Minier, V., et al. 2014, A&A, 564, A106

Urquhart, J. S., Figura, C. C., Moore, T. J. T., et al. 2015, MNRAS, 452, 4029

Urquhart, J. S., Morgan, L. K., Figura, C. C., et al. 2011, MNRAS, 418, 1689

Uyaniker, B., Fürst, E., Reich, W., Aschenbach, B., & Wielebinski, R. 2001, A&A, 371, 675

Vázquez-Semadeni, E., González-Samaniego, A., & Colín, P. 2017, MNRAS, 467, 1313

Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, submitted (arXiv:1903.11247)

Vázquez-Semadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245

Walsh, J. R., Ogura, K., & Reipurth, B. 1992, MNRAS, 257, 110

Wang, H., Yang, J., Wang, M., & Yan, J. 2003, AJ, 125, 842

Wang, J., Townsley, L. K., Feigelson, E. D., et al. 2008, ApJ, 675, 464

Ward-Thompson, D., Kirk, J. M., André, P., et al. 2010, A&A, 518, L92

Werner, M. W., Becklin, E. E., Gatley, I., et al. 1979, MNRAS, 188, 463

White, G. J., Nelson, R. P., Holland, W. S., et al. 1999, A&A, 342, 233

Wienen, M., Wyrowski, F., Schuller, F., et al. 2012, A&A, 544, A146

Williams, R. J. R., Ward-Thompson, D., & Whitworth, A. P. 2001, MNRAS, 327, 788

Willson, R. F., & Folch-Pi, F. J. 1981, AJ, 86, 1084

Wilson, B. A., Dame, T. M., Masheder, M. R. W., & Thaddeus, P. 2005, A&A, 430, 523

Wilson, T. L., & Mauersberger, R. 1990, A&A, 239, 305

Wolf-Chase, G. A., Walker, C. K., & Lada, C. J. 1995, ApJ, 442, 197

Wright, N. J., Parker, R. J., Goodwin, S. P., & Drake, J. J. 2014, MNRAS, 438, 639

Wu, J., Evans, N. J., II, Gao, Y., et al. 2005, ApJL, 635, L173

Wu, J., Evans, N. J., II, Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313

Xu, Y., Moscadelli, L., Reid, M. J., et al. 2011, ApJ, 733, 25