The Bolocam Galactic Plane Survey. XIV. Physical Properties of Massive Starless and Star Forming Clumps

We sort $4683$ molecular clouds between $10^\circ<\ell<65^\circ$ from the Bolocam Galactic Plane Survey based on observational diagnostics of star formation activity: compact $70$ $\mu{\rm m}$ sources, mid-IR color-selected YSOs, ${\rm H_2O}$ and ${\rm CH_3OH}$ masers, and UCHII regions. We also present a combined ${\rm NH_3}$-derived gas kinetic temperature and ${\rm H_2O}$ maser catalog for $1788$ clumps from our own GBT 100m observations and from the literature. We identify a subsample of $2223$ ($47.5\%$) starless clump candidates, the largest and most robust sample identified from a blind survey to date. Distributions of flux density, flux concentration, solid angle, kinetic temperature, column density, radius, and mass show strong ($>1$ dex) progressions when sorted by star formation indicator. The median starless clump candidate is marginally sub-virial ($\alpha \sim 0.7$) with $>75\%$ of clumps with known distance being gravitationally bound ($\alpha<2$). These samples show a statistically significant increase in the median clump mass of $\Delta M \sim 170-370$ M$_\odot$ from the starless candidates to clumps associated with protostars. This trend could be due to (i) mass growth of the clumps at $\dot{M}\sim200-440$ Msun Myr$^{-1}$ for an average free-fall $0.8$ Myr time-scale, (ii) a systematic factor of two increase in dust opacity from starless to protostellar phases, (iii) and/or a variation in the ratio of starless to protostellar clump lifetime that scales as $\sim M^{-0.4}$. By comparing to the observed number of ${\rm CH_3OH}$ maser containing clumps we estimate the phase-lifetime of massive ($M>10^3$ M$_\odot$) starless clumps to be $0.37 \pm 0.08 \ {\rm Myr} \ (M/10^3 \ {\rm M}_\odot)^{-1}$; the majority ($M<450$ M$_\odot$) have phase-lifetimes longer than their average free-fall time.


INTRODUCTION
Massive stars (M > 8 M ) strongly influence the evolution of galaxies and the ISM, yet their formation remains an open problem in contemporary astrophysics (see reviews by McKee & Ostriker 2007;Tan et al. 2014). A major bottleneck in the study of massive star formation is the difficulty in identifying and systematically analyzing the properties of the incipient phases of intermediate-and high-mass (protocluster) star formation, high-mass starless clumps. These objects fragment into dense starless cores that subsequently contract to form individual or bound systems of protostars . Starless cores (and gravitationally bound prestellar cores) are cold (T d < 20 K), dense (n(H 2 ) > 10 4 cm −3 ), and centrally concentrated objects that lack an embedded protostar (di Francesco et al. 2007). The study of low-mass (< 10 M ) starless cores has received considerable attention due to their proximity and the systematic mapping of nearby molecular clouds in the Gould's Belt (Ward-Thompson et al. 1994;Tafalla et al. 2004;André et al. 2007André et al. , 2013. In contrast, the study of intermediate and high-mass starless cores and clumps is still in a nascent state due in part to their greater distances and few systematic observations to identify them (Wilcock et al. 2012;Tackenberg et al. 2012;Ragan et al. 2013;Traficante et al. 2015a). Our understanding of how cluster formation is initiated and ensuing protocluster evolution both ultimately depend on identifying and constraining the physical properties of representative samples of high-mass starless clumps. This situation is beginning to change with the publication of a catalogs of 170 and 667 starless clumps identified with the Herschel Space Observatory towards Infrared Dark Clouds (IRDC;Wilcock et al. 2012;Traficante et al. 2015a). In this paper, we analyze the star formation activity of clumps observed in the Bolocam Galactic Plane Survey to identify a population of over 2000 massive starless clump candidates in a 84.4 square degree region of the first quadrant of the Galaxy.
Recent blind surveys of the dust continuum emission at (sub-)millimeter wavelengths are able to detect embedded star forming regions in a wide range of evolutionary states, including the starless phases (e.g., Tackenberg et al. 2012;Traficante et al. 2015a), throughout the Galaxy. There are three primary surveys that have mapped far-infrared through millimeter wavelength emission in the Milky Way. The Bolocam Galactic Plane Survey 10 (BGPS; Aguirre et al. 2011;Rosolowsky et al. 2010;Ginsburg et al. 2013) at λ = 1.1 mm surveyed from −10 • < < 90 • and targeted regions in the second and third quadrants with |b| < 0.5 • (expands to |b| < 1.5 at selected ) at 33 resolution. ATLASGAL (Schuller et al. 2009;Contreras et al. 2013;Csengeri et al. 2014) surveyed at λ = 870 µm between −60 • < < 60 • and |b| < 1.5 • with an extension covering −80 • < < −60 • and −2.0 < b < +1.0 at 22 resolution. In the far-infrared and submillimeter, the Herschel Infrared Galactic Plane Survey (Hi-GAL; Molinari et al. 2010) surveyed the entire Galactic plane in the five PACS and SPIRE bands from λ = 60 − 500 µm. In this paper, we shall focus on studying the physical properties of starless and star forming clumps identified in the BGPS.
The BGPS survey identified 8294 clumps with 1.1 mm flux densities above 100 mJy and 8594 clumps cataloged in total (Rosolowsky et al. 2010;Aguirre et al. 2011;Ginsburg et al. 2013). The second data release (v2.0) of the Bolocam images and source catalog (Ginsburg et al. 2013) used improved analysis of the time-stream data to improve angular flux recovery with full recovery out to ∼ 80 and partial recovery out to ∼ 300 . The v2.0 data also resolves systematic flux and pointing calibration offsets and includes new fields in the second quadrant. In this work we exclusively use the BGPS v2.0 data products, and all clumps are referenced to the BGPS v2.0 maps and Bolocat clumps.
Extensive followup observations of BGPS clumps have been carried out to determine unique velocities, kinematic distances, and gas kinetic temperatures to use in calculations of physical properties. Followup spectroscopic observations of the BGPS catalog in dense molecular gas tracers HCO + 3 − 2 and N 2 H + 3 − 2 with the Arizona Radio Observatory's Heinrich Hertz Submillimeter Telescope (HHSMT; Schlingman et al. 2011;Shirley et al. 2013) detected about half of the objects in the BGPS catalog (> 97% of those clumps have a single velocity component) from which kinematic distances may be derived. A novel Bayesian technique has been developed to calculate a continuous function of the probability that a source lies at a given heliocentric distance (Ellsworth-Bowers et al. 2013. The posterior Distance Probability Density Function (DPDF) is constructed by multiplying the bimodal likelihood function (with equal probability of being located at the near or far kinematic distance) with prior distributions that account for the source positions relative to the known H 2 distribution in the Galaxy and the coincidence and contrast of 8 µm absorption features (Ellsworth-Bowers et al. 2013). Since the dense molecular gas lines are only detected toward approximately half of BGPS clumps, Ellsworth-Bowers et al. (2014) extended the number of sources with bimodal likelihood functions by developing a morphological matching technique between 13 CO 1 − 0 from the Galactic Ring Survey (Jackson et al. 2006) and the BGPS 1.1 mm emission. This technique assigns a clump the dominant 13 CO velocity component that matches the 1.1 mm morphology. The DPDF formalism is a powerful tool that permits proper accounting and propagation of the distance uncertainty in calculations of distance dependent physical properties (i.e., size, mass, and luminosity).
Early BGPS papers have characterized clump physical properties on modest subsets of the BGPS. Schlingman et al. (2011) determined the physical properties of 529 sources drawn uniformly from the BGPS 1.1 mm flux distribution, finding that BGPS sources have median physical radii of 0.75 pc, median total masses of 330 M , and non-thermal linewidths that are 10 times the thermal linewidth on median. Dunham et al. (2011b) performed an initial evolutionary analysis on a subsample of 456 BGPS sources towards which NH 3 was detected, finding a median gas kinetic temperature of 16 K and identifying a subset of starless clump candidates (SCCs) that composes ∼ 1/3 of their sample. While BGPS sources are commonly referred to generically as clumps, Dunham et al. (2011b)  Observed S1.1, ¢v, µ eq , ::: Physical Properties : M, §, R eq , ®, ::: Statistical Quantities : ¹1/2, ¹, ¾, ::: which sample cores, clumps, and clouds on different scales and at different distances. Indeed, higher resolution observations at 350 µm indicate that BGPS clumps typically break up into collections of smaller substructures (Merello et al. 2015). All of these early studies made crude assumptions to resolve the kinematic distance ambiguity (i.e., all sources associated with 8 µm absorption are at the near distance or using the presence or lack of HI self-absorption to distinguish between near and far kinematic distance; see §4.3). Ellsworth- Bowers et al. (2015) is the first study to utilize the full power of the DPDF formalism in which they characterize the clump mass distribution for the subset of 1710 clumps with well-constrained DPDFs. They find that the clump mass distribution is best fit with a lognormal distribution with an approximately powerlaw distribution for high-mass sources and a power-law index that is intermediate between that observed for giant molecular clouds and that observed for the stellar initial mass function.
A significant remaining challenge with the BGPS data set is to identify the star formation activity of all of the BGPS clumps and, in particular, identify the earliest phase of protocluster formation: starless clumps. The extensive characterization and followup of the BGPS presents an excellent sample drawn from a blind survey to rigorously calculate population statistics and physical properties. Measuring the physical properties of starless clumps will not only constrain the initial conditions of protocluster evolution but starless clumps will also be the proving grounds between theories of high-mass star formation, such as virialized turbulent core accretion (TCA; McKee & Tan 2003) and a sub-virial competitive accretion (CA; Bonnell et al. 2001Bonnell et al. , 2007. Resolving the current tension between these two theories is a major objective in the near-term study of high-mass star formation, but requires the identification of a representative sample of massive starless clumps for further high resolution study. In this paper we address questions related to how starless clumps, as the host environments of future intermediate-and high-mass stars, evolve through the starless phase into the protostellar phase. How long does the starless phase last and how does the starless clump lifetime depend on the the mass of the clump? Do clumps undergo significant evolution in their fundamental physical properties and do they undergo significant interaction with the surrounding molecular cloud throughout their lifetimes? We study these questions by analyzing large, statistical subsamples (N ∼ 10 2 − 10 3 ) of BGPS clumps.
We present a catalog of 4683 BGPS clumps sorted by observational indicators of star formation activity from a comparison to other Galactic plane surveys in a common overlap range between 10 • < < 65 • . We also present 1215 targeted NH 3 and H 2 O maser observations with the 100 m Robert C. Byrd Green Bank Telescope (GBT) in §2. Schematically, we follow the procedure shown in Figure 1. In §3 we describe the survey data sets and methods in developing a BGPS star formation catalog. In §4 we describe the methods used to extend the sample of Distance Probability Density Functions to additional sources with a position-position-velocity (PPV) clustering algorithm. We use Monte Carlo sampling of the DPDFs, re-sampling techniques to account for contamination of evolved stars, and sampling of uncertainties on observed quantities to calculate and analyze the physical properties of the BGPS clumps associated to star formation indicators in §5. We discuss the properties of the starless clump candidate sample and possible evidence for clump mass growth during this phase in §6.

NH 3 AND H 2 O MASER SURVEY OF BGPS SOURCES
2.1. New GBT NH 3 Observations We have conducted targeted, spectroscopic observations of the NH 3 (1, 1), (2, 2), and (3, 3) inversion transitions as well as the J K+,K− = 6 1,6 − 5 2,3 22 GHz H 2 O maser line with the GBT 11 toward 1215 new BGPS sources. Dunham et al. (2011b) observed 631 BGPS sources within narrow ranges in Galactic longitude, while the new data presented in this work observed BGPS clumps that were detected in dense gas tracer HCO + 3 − 2 and had T pk (HCO + ) > 0.3 K from Schlingman et al. (2011) and Shirley et al. (2013). All new observations were targeted with the center pixel of the K-band Focal Plane Array (KFPA) towards the S 1.1mm flux density peak in the BGPS v2.0 maps (Ginsburg et al. 2013). We measure NH 3 properties using the slab model previously adopted in Dunham et al. (2011b). This model assumes a uniform-temperature, beam-filling slab of NH 3 with column density N (NH 3 ), intrinsic velocity dispersion σ v , kinetic temperature T K , and line-of-sight velocity v LSR . The model assumes equal column densities of ortho-and para-NH 3 and that the levels are populated in thermodynamic equilibrium. Given these assumptions, we use a non-linear least square regression to determine the optimal emission model incorporating the hyperfine structure 11 The GBT is operated by the The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. and opacity of each component that matches the observed data (Rosolowsky et al. 2008). As in previous pointed spectroscopy studies, the NH 3 fitting provides an excellent reproduction of the observed data with the derived parameters representing emission-weighted averages of the physical conditions within the GBT beam. The observational setup, data reduction methods, and NH 3 fitting techniques are described in detail in Dunham et al. (2010Dunham et al. ( , 2011b. Previous KFPA observations based on the BGPS v1.0 are updated to match the v2.0 based on the coincidence of the pointing within the BGPS label-mask (discussed in §3.1 describing the catalog cross matching method). Because of the changes in the BGPS catalog between data releases, 178 pointings are duplicates within the same v2.0 clumps and 245 pointings are not within a half-beam of a v2.0 clump. We use the observation closest to the peak flux density S 1.1mm position when there are multiple NH 3 observations within the clump. Columns are included in Table 2 for the Bolocat v2.0 source and coordinate offset from the peak flux position. Table 2 also includes the observations from Dunham et al. (2011b) for consistency with the transition to BGPS v2.0. We find an overall detection rate in NH 3 (1,1) at 75%. For accurate estimates of the kinetic temperature we require a source to have > 3σ detections in both the (1,1) and (2,2) transitions, resulting in 1663 clumps with suitable T K fits. Table 2 lists the observed NH 3 properties from the fits and 3σ upper limits for non-detections. Unique T K measurements are added for 48 BGPS clumps in 10 • < < 65 • from observations in Wienen et al. (2012).

New GBT H 2 O Maser Observations
The 22 GHz H 2 O maser line is a well known tracer of protostellar activity and was simultaneously observed in the fourth IF alongside the NH 3 (1, 1) − (3, 3) inversion transitions. A separate reduction in GBTIDL is used for the H 2 O maser data because of the complex structure in the spectra from multiple overlapping maser spots in the beam, strong baseline ripples, and a telluric feature at the topocentric zero velocity. The observations are dual-frequency switched with a 5 MHz (69 km s −1 ) throw. The telluric feature is pressure broadened by ∼ 2 MHz (∼ 30 km s −1 ) and adds significant structure to the baseline when frequency switched. We baseline fit the spectra with first to third order polynomials to accommodate the complex baseline morphology. Spectra are required to have a > 5σ detection in a single channel to qualify as detections. Table 3 lists the observed H 2 O maser properties and 3σ upper limits for non-detections. These values include the R.A. and Dec. of the GBT pointing, offset from the BGPS peak position θ offset , velocity of the peak component v LSR , FWZI spread between the left-and right-most velocity components v spread , mainbeam peak intensity T mb , main-beam integrated intensity W , and number of unblended velocity components N lines . We observe 472 detections, 343 of which are unique associations not contained in the published surveys by the Mopra and Arcetri observatories discussed in §3.3.1. No significant deviation was seen between the center velocity defined by the first moment and that of the strongest peak. The number of unblended line peaks are identified through visual inspection, but because the maser components can be heavily blended towards strong sources, this number represents a lower limit to the number of maser components towards a BGPS source. For extreme H 2 O maser emission, the weaker telluric feature is covered and adds a ∼ 0.5 K km s −1 uncertainty to the calculated H 2 O maser integrated intensity. This overlap introduces a < 1% uncertainty in the > 100 K km s −1 H 2 O maser spectra where this occurs.

DEVELOPING A STAR FORMATION INDICATOR
CATALOG Recently completed blind Galactic plane surveys with wavelength coverage from the near infrared to radio provide a wide range of indicators of star formation activity. We develop an evolutionary catalog for clumps in the BGPS based on the indicators of star formation that can be associated with the clumps. These evolutionary indicators are then grouped into the categories so that their observed and physical properties can be compared. From this catalog, we identify a subsample of starless clump candidates that lack indicators of protostellar activity. Below we describe the surveys and data sets used to define the flags of star formation indicators. The published surveys include: the Red MSX Survey (Lumsden et al. 2002), Glimpse Red Source Catalog (Robitaille et al. 2008), Hi-GAL 70 µm images (Molinari et al. 2010), H 2 O Southern Galactic Plane Survey (Walsh et al. 2011), Methanol Multibeam Maser Survey (Breen et al. 2015, and Coordinated Radio and Infrared Sky Survey for High-mass Star Formation (Hoare et al. 2012) in addition to the GBT H 2 O maser and NH 3 observations presented in this paper. We aim for the most uniform coverage in the surveys to accurately compare star formation indicator populations. The common overlap region among the surveys described below is between 10 • < < 65 • .

Catalog Cross-Matching
All published catalogs are cross-matched to the BGPS v2.0 catalog based on line-of-sight coincidence on the sky. The Bolocat modified seeded watershed algorithm extracts the clump catalog from the Bolocam data and defines the spatial extent of each clump (Rosolowsky et al. 2010). In this context the "label-map" or "label-mask" refers to the pixels in the Bolocam 1.1 mm flux image associated with a clump by the watershed. For a source from another catalog to be matched with a clump we require that it must be contained within the clump label-mask. In the common overlap region between 10 • < < 65 • there are 4683 BGPS clumps and 2203 clumps with either HCO + or NH 3 dense gas detections. Table 1 lists the observed 1.1 mm properties for the 4683 BGPS clumps within the common overlap region. Table 4 shows the detection statistics for the cross-matched catalogs described in §3.2.1- §3.4 to the BGPS. The table lists the number of sources in the overlap region, the number of sources matched with BGPS clumps, and the number of BGPS clumps containing at least one matched source. Without accurate independent distances for both a BGPS clump and a catalog source for a star formation tracer, we must rely on spatial coincidence on the sky to cross-match sources. To assess the accuracy of the cross-matched associations, Dunham et al. (2011a) estimate the rate of chance alignments statistically. Dunham et al. (2011a) compare the fraction of cross-matched sources in a given catalog to the average fraction of sources cross-matched when the coordinates have been randomized on area much greater than a single BGPS clump's environment (see Table 1). Randomizing the coordinates of the catalog preserves the large-scale projected spatial distribution of sources throughout the Galaxy. If the sources in a catalog are quasi-uniformly distributed on the sky, this randomized overlap fraction would simply be the area of the sky covered by BGPS clumps divided by the survey coverage. However, because of Galactic structure and line of sight projection effects, the distribution of most catalog sources (e.g., red GLIMPSE sources) will not be uniform, and this will be reflected in a higher count of randomized cross-matches towards, e.g., the inner Galaxy. The randomized cross-matching yields typically between 7 − 9% of catalog sources being associated to BGPS clumps. In this work we use the BGPS v2.0 whereas Dunham et al. (2011a) used the BGPS v1.0. We find similar results for the RMS, EGO, and R08 sources matched to clumps ( §3.2) and the number of clumps with associated sources within a few percent of the total fraction of matched sources. Catalogs that are new additions in this work such as MMB and COR-NISH also have high-association rates, > 80%, so we do not suspect that the sources in these catalogs are highly contaminated by chance alignments. Furthermore, other datasets in this work such as the Hi-GAL visual inspection and the GBT H 2 O observations are targeted and the coordinates can not be randomized.

Red MSX Survey
The Red MSX Survey (RMS; Lumsden et al. 2002Lumsden et al. , 2013 with observations at 8, 12, 14, and 21 µm defines a sample of massive YSOs complete above L bol > 10 4 L and d < 10 kpc, approximately a B0 type main sequence star. The RMS catalog region has complete overlap with the BGPS for > 10 • , where confusion does not dominate the MSX data. Extensive followup observations of RMS sources yield detailed classifications: proto-planetary nebula, planetary nebula, evolved star, HII region, YSO, OH/IR star, young/old star, HII/YSO, carbon star, and a miscellaneous "other". We use the RMS classifications associated with YSOs (YSO, HII/YSO, young/old star) to classify BGPS sources in the mid-IR category. We do not associate BGPS sources with star formation activity for the other categories because contamination between the RMS classifications are low. This is further evidenced by the minimal overlap in Table 4. The RMS sources provide a well-vetted sample of massive YSOs, but for spectral types later than B0, or YSOs that are at distances > 10 kpc, the survey is incomplete. Recent Galactic plane surveys by Spitzer and Herschel offer sensitive datasets through the midand far-IR to probe the star formation activity in BGPS clumps and complement the RMS data.

GLIMPSE Red Source Catalog
The GLIMPSE survey is a Spitzer/IRAC Legacy survey of the Galactic mid-plane, at 3.6, 4.5, 5.8, 8.0 µm (Benjamin et al. 2003;Churchwell et al. 2009). The GLIMPSE I and II fields provide complete overlap coverage with the BGPS for 10 • < < 65 • . Robitaille et al. (2008) (hereafter R08) use the GLIMPSE Catalog to a select a total number of 18,949 intrinsically red ([4.5] − [8.0] 1) sources that also meet photometric quality criteria. MIPS 24 µm magnitudes are also calculated at the same positions for 86.9% of the sources from the MIPS-GAL survey (Carey et al. 2009). Because GLIMPSE has better spatial resolution and sensitivity than the MSX-survey, lower luminosity YSOs can be detected. R08 estimate that a 100 L Stage I YSO should be detected at distances between 0.2 − 10 kpc and a 10 4 L YSO between 2 − 30 kpc. Intrinsically red sources are further classified by color-magnitude and colorcolor selection criteria (see eq. 7-9 Robitaille et al. 2008) into candidate YSOs, "Standard" C-and O-rich AGB stars (sAGB), and "Extreme" AGB stars (xAGB). The color criteria selects sources with spectral index α −1.2, which following the "Class" system of Greene et al. (1994) selects all Class I and a significant number of Class II YSOs. Contamination from planetary nebula and background galaxies is expected to be between 2 − 3%. However, because AGB and YSO sources have similar mid-IR colors, the YSO/AGB discrimination is uncertain based on 4.5, 8.0, and 24 µm photometry alone, with up to 50% contamination between the two categories. Dunham et al. (2011b) estimate a ∼ 40% contamination fraction of YSOs in the R08 candidate AGB sources. We shall account for this contamination in the Monte Carlo sampling in §5. The number of R08 candidate YSO or AGB sources matched to BGPS clumps is listed in Table 4. An additional class of sources called Extended Green Objects (EGOs) have been identified from the Spitzer GLIMPSE data (Cyganowski et al. 2008;Chen et al. 2013). EGOs are sources of extended 4.5 µm emission from shocked H 2 driven by strong outflows from massive star formation. 72 BGPS clumps match EGOs, and all except one clump (G044.008-00.026) contain RMS YSOs, R08 YSOs, or 70 µm compact sources (see §3.2.4).

MIPSGAL 24 µm
The Version 3.0 data products from the Spitzer MIPSGAL survey cover −69 • < < 68 • and |b| < 1 • at 24 µm (Carey et al. 2009). The MIPS 24 µm band is sensitive to the hot dust found around YSOs and HII regions, but also to envelopes of evolved stars. Gutermuth & Heyer (2015) have created a comprehensive point-source catalog from the MIPSGAL 24 µm maps. Two catalogs are presented: a high-quality catalog and more complete archive. Combined the catalogs total 268,430 point sources in the fields of the BGPS overlap region, with 4.9% of MIPSGAL sources associated to clumps and 80.9% of clumps associated to MIPSGAL point sources. While some subset of MIPSGAL sources do represent observational indicators of embedded protostellar activity, this low match fraction of 4.9% is consistent with the 6.8% ± 2.2% randomized match fraction calculated in Dunham et al. (2011b) for RMS Evolved Stars (i.e., sources with no physical association). Figure 2 illustrates this poor correspondence to BGPS clumps. Two maps at 70 µm and 24 µm are shown of the same field centered on a single, isolated clump. Whereas the 24 µm image shows identified point sources from both the Catalog and Archive, only a few point sources have direct association to the clump 1.1 mm emission. In contrast, the 70 µm image uniquely associates a single compact source to the clump. The 70 µm dust continuum emission is more sensitive to the deeply embedded star formation activity that is better positionally correlated to BGPS clumps. We conclude that line of sight association of 24 µm point sources from the MIPSGAL catalog are not a sufficient sole indicator of star formation activity because of the high contamination from evolved stars. Thus while the MIPSGAL flags are not used in any calculations of physical properties, for completeness we present the flags in Table 4 (see §3.5). These contamination issues also strongly affect the lower resolution WISE 22 µm observations. To further refine our set of star formation indicators, rather than pursue colorand quality-selection criteria to identify a less contaminated subsample from the MIPSGAL catalog, we associate clumps to 70 µm compact sources.

Herschel Hi-GAL 70µm
The Herschel Infrared Galactic Plane Survey (Hi-GAL) is a far-IR survey of the Galaxy with the 70 and 170 µm PACS bands and the 250, 350, and 500 µm SPIRE bands (Molinari et al. 2010). The Hi-GAL survey completely covers the BGPS fields that we discuss. Previous survey data sets in the 100 µm regime from IRAS and Akari are not suitable for detecting faint point sources due to their poor angular resolution (> 1 ). The Herschel Hi-GAL 10.2 resolution (Traficante et al. 2011) at 70 µm is essential to disentangling the far-IR emission reprocessed by dust around embedded YSOs and their often active environment. Most important however is the less severe contamination in the Hi-GAL 70 µm images from evolved stars compared to the MIPS 24 µm or WISE 22µm images. Whereas > 80% clumps have at least one associated MIPSGAL source, only 46% of clumps have an associated 70 µm compact source from visual identification with clearly correlated positions to the 1.1 mm maps.
Due to complicated background structure, identifying compact sources algorithmically is difficult. The CuTeX algorithm (Molinari et al. 2011) was developed specifically for identifying sources in highly spatially variable backgrounds such as those found in the Hi-GAL images. Careful visual inspection is necessary because of false negatives where CuTeX misses faint sources and false positives where CuTeX mis-identifies sources in regions of strong emission. Cuts on CuTeX fit parameters were trained on publicly available data in the = 30 • field  and then applied to all cutout images to compare with the results of the visual inspection classifications. The comparisons show that 77.7% of the classifications agree, and 22.3% disagree. 81% of the disagreements originate from sources flagged with fainter 70 µm sources from the visual inspection but missed by CuTeX.
We visually inspect 5 × 5 cutouts of the reduced Hi-GAL 70 µm band images (Molinari et al. 2011) for comparison against the BGPS data. For consistency with the catalog matching scheme used for the other surveys in this study, we classify clumps as containing a Hi-GAL 70 µm source if it is within the Bolocat label-mask. Visual inspection was carried out only on the 70 µm images alone, without simultaneous comparison to GLIMPSE or MIPSGAL images to prevent systematic bias in the classifications. We assign BGPS clumps one of five flags in the visual inspection process. Multiple flags are necessary to qualify the image background level and source characteristics. Flags 1, 2, and 4 represent clumps identified with a compact source within the label-mask while flags 0 and 3 lack a clearly identifiable source. Representative 70 µm cutouts of these flags are shown in Figure 3. Flag 1's show a high confidence compact source on the scale of the Herschel 10.2 PSF. Flag 2's are more diffuse than the PSF but identifiable as compact sources due to strong, localized emission above the background level. These are common for extremely bright 70 µm sources with 80% of uniquely high-mass indicators, RMS sources, CH 3 OH masers, and UCHII regions, assigned flag 2's. We identify 944 flag 1 and 1044 flag 2 clumps, or 42.5% of clumps. Flag 4's show a weak, PSF-like source whose classification is lower confidence. Adding the 182 flag 4's increases the total to 46.3% of clumps classified with a Hi-GAL source. The brightness and complexity of the 70 µm background can make the visual source identification difficult. Because of this we assign 1847 flag 0's to clumps with quiescent, smooth, and low-level backgrounds and 655 flag 3's to clumps with bright and more complex backgrounds. The far-IR background of the Milky Way changes locally and as a function of Galactic longitude (generally 10 2 − 10 3 MJy sr −1 ), and thus the flux threshold for source detection changes as well.  The high rate of overlap between 70 µm sources and more extreme indicators along with the large number of unique 70 µm associations suggest that the Hi-GAL data is broadly probing a lower YSO bolometric luminosity, L bol , than the other indicators. To assess the low-luminosity limit of the associations, we compute L bol completeness functions for a representative sample of background levels given the sample of DPDFs in the Distance Catalog. The 70 µm background level and complexity is the limiting factor in the identification of compact sources. A detailed analysis of the background level distribution for the full Hi-GAL 70 µm dataset with the identification threshold quantified through fake source injection is beyond the scope of this paper. Instead, we bracket the completeness function with a sample of 30 clumps associated to the least confident 70 µm identifications (classified as Hi-GAL flag 4), evenly split between low background levels of 500 − 1000 MJy sr −1 and high background levels of 2000 − 3000 MJy sr −1 . Values are calculated from the cutouts described above. Aperture photometry of the 70 µm sources in these selected regions yields sources identified with flux densities 0.3 Jy for the low backgrounds and 1 Jy for the high backgrounds. Using these two detection thresholds, we compare to fluxes derived from the 20,000 axisymmetric YSO radiative transfer models computed in Robitaille et al. (2006Robitaille et al. ( , 2007. We use the PACS 70 µm fluxes from the largest physical aperture, 100,000 AU, which corresponds to the 20 diameter aperture at the median clump distance of 5 kpc (see §5.1.2). We then draw from Monte Carlo simulations for each clump distance probability density function (described in §5.1) to scale the Robitaille et al. (2006) PACS 70 µm model flux and apply the detection threshold cut. Figure 4 shows the completeness function for all clumps with well-constrained distances and the cyan line shows the same but for the subset of clumps between 28 • < < 31 • for which publicly available Hi-GAL data exists ). The dotted lines show the 50% and 90% limits. We find in the low background environments the 70 µm visual flags are 50% complete to YSOs with L bol > 10 L and 90% complete to YSOs with L bol > 50 L . For the high background environments, we find the 50% completeness limit to YSOs at L bol > 45 L and 90% complete to YSOs with L bol > 140 L . To put these luminosities in perspective, a ∼ 1 M and ∼ 5 R YSO accreting at ∼ 10 −5 M yr −1 corresponds an accretion luminosity of ∼ 30 L (N.B. episodic accretion rates result in large variations in the accretion luminosity at a given epoch; see Dunham et al. 2014). While the true sample completeness limits likely lie between these two values, it should be kept in mind that these limits are below the L bol of the other star formation indicators by several orders of magnitude.
Among the positive Hi-GAL classifications, 1027/2170 (47.3%) have one or more additional star formation indicator and 1143/2170 (52.7%) are Hi-GAL 70 µm unique, representative a sample of deeply embedded and early protostellar candidates. In the following analysis, this sample of clumps uniquely identified with star formation activity from 70 µm compact sources is treated separately as its own category. It should be noted that the 70 µm unique category sources may have a MIPS 24 µm counterpart but does not contain any other star formation indicator (unique is defined in terms of our set of star formation indicators).

H2O Masers
We complement the mid-and far-IR star formation indicators with radio observations of masers. The outflows driven by young protostars shock dense regions of the ISM and collisionally pump H 2 O that creates a population inversion which drives the maser emission. The typical physical conditions for a masing regions are a spatial density n ∼ 10 7 − 10 9 cm −3 , kinetic temperature T K ∼ 400 K, spatial extent of d ∼ 0.66 AU, and an aspect ratio ∼ 50 : 1 (Hollenbach et al. 2013). H 2 O masers are time variable on scales of months to years. The beaming effect of H 2 O masers make them a sensitive probe at large distances but because of the necessary narrow geometrical alignment, the non-detection of a H 2 O maser does not imply the lack of star formation activity. In addition to YSO molecular outflows, H 2 O masers are also observed in the shells around AGB stars (Imai et al. 2002). From high resolution followup of the H 2 O Southern Galactic Plane Survey (HOPS; Walsh et al. 2011) and association to catalogs of star formation indicators Walsh et al. (2014) find that 69% of water maser sites are associated with star formation, 19% with evolved stars, and 12% with unknown identification. We expect limited contamination of AGB stars coincident with our GBT pointings because the observations are targeted towards the clump peak S 1.1mm and N (H 2 ) positions (θ HPBW ≈ 30 ) where a YSO is likely to be located, but a chance alignment by a field AGB star is unlikely. Thus we use the detection of a H 2 O maser as an observational indicator of star formation activity.
Large observational H 2 O maser survey programs have been  carried out by the Arcetri survey (Valdettaro et al. 2001), HOPS (Walsh et al. 2011), and the RMS H 2 O Maser Survey ). The variability of H 2 O masers make multiple catalogs at different epochs complementary for systems that drop below a surveys sensitivity limits (Furuya et al. 2003). The Arcetri catalog aggregates observations of all known literature H 2 O water masers prior to 2001. HOPS is a blind Galactic plane survey with overlap coverage with the BGPS extending to < 20 • . The HOPS team discovered a number of new H 2 O masers not found in the Arcetri catalog, most in the southern sky that is inaccessible to the Arcetri survey. While the survey overlap between HOPS and the BGPS-GBT H 2 O maser survey is limited to 10 • < < 20 • , it provides valuable context to evaluate the targeted GBT observations. In this range, the GBT observations have a 3.4 times higher detection rate towards BGPS clumps. This is likely due to the tighter beam coupling of the GBT compared to Mopra and the lower baseline RMS for the targeted, single-pointing observing strategy. The GBT observations primarily detect weaker H 2 O masers with median integrated intensity ∼ 5 Jy km s −1 compared to HOPS with ∼ 5 × 10 2 Jy km s −1 .

CH3OH Masers
The 6.7 GHz Class II methanol (CH 3 OH) maser is thought to be exclusively associated with high-mass star formation because of the environments necessary for IR pumping (Breen et al. 2013). A comparison between 6.7 GHz CH 3 OH maser emission from the MMB and ATLASGAL 870 µm dust emission finds that 99% of MMB sources are associated with AT-LASGAL emission indicating that methanol maser emission is almost ubiquitously associated with massive protostellar clumps (Urquhart et al. 2015). We associate BGPS clumps to CH 3 OH maser indicators based on available literature data. We use data from the Arecibo CH 3 OH survey (Pandian et al. 2007(Pandian et al. , 2011, the Methanol Multibeam Survey (MMB; Breen et al. 2012bBreen et al. ,a, 2014Breen et al. , 2015, and aggregated literature data in Pestalozzi et al. (2005). Within the overlap region, MMB data extends from 10 • < < 60 • , Arecibo data extends from 35 • < < 54 • , and the data in Pestalozzi et al. (2005) are all-sky targeted observations. Thus there is a small gap in the CH 3 OH maser survey data from 60 • < < 65 • , but there are only 15 BGPS clumps in this range, so there should not exist a strong sampling bias. While the MMB is the most uniform and complete of the surveys, multiple overlapping datasets are useful because the intensity of CH 3 OH masers is time variable. This can be important for masers that drop below the sensitivity limits of a single survey. With multiple surveys at different epochs we create a more complete catalog of CH 3 OH masers.

Radio Continuum -UCHII Regions
UCHII regions are an unambiguous indicator of high-mass star formation because an OB star must be present to internally photoionize the clump. We compare the BGPS with the Coordinated Radio and Infrared Sky Survey for High-Mass Star Formation (CORNISH; Hoare et al. 2012;Purcell et al. 2013) observed with the JVLA in 5 GHz continuum from 10 • < < 65 • . We compare all 241 sources classified as UCHII region candidates in Purcell et al. (2013), excluding 43 sources classified as either "HII-Region", "Diffuse HII-Region", or "Dark HII-Region", representative of later stages of evolution. We find 219/243 (90.1%) CORNISH UCHII candidates associated with 170 BGPS clumps. No BGPS clumps are exclusively identified with star formation activity from the CORNISH UCHII region sources. All such clumps are associated with an additional evolutionary indicator: 169 Hi-GAL 70 µm, 69 mid-IR, 103 H 2 O maser, and 70 CH 3 OH maser. Urquhart et al. (2013b) compare CORNISH with ATLASGAL and find 213/243 are matched to 170 ATLASGAL clumps. The 24 CORNISH sources classified as UCHII region candidates that are not associated with millimeter continuum emission in the BGPS could potentially be BGPS non-detections, contamination/mis-identification in the CORNISH survey, or more evolved systems that have disrupted most of the surrounding clump.

Star Formation Indicator Statistics
We combine these catalogs of observational indicators of star formation activity into a series of flags to sort clumps into groups for joint analysis of their physical properties. The most basic distinction is between clumps that are associated with no star formation indicators, the "Starless Clump Candidates" (SSC), and "Protostellar" Clumps that are associated to one or more indicator. The star formation indicator groups are then ranked by the approximate luminosity of a single protostar that produces the indicator (see Figure 5). These rankings are a heuristic motivated by proposed linear clump evolutionary sequences (Battersby et al. 2010(Battersby et al. , 2011Chambers et al. 2009), but the ordering is settled upon because it consistently yields monotonic trends for clump properties. "70 µm Unique" are clumps that only contain a Hi-GAL compact source, and thus are candidates for deeply embedded (≈ Class 0/I) protostellar activity. "Mid-IR" clumps are associated to one or more R08 YSO-flagged source, RMS protostellar source, and/or EGO. We group R08, RMS, and EGOs together because they are frequently coincident with the same protostellar sources and they use similar wavelengths to identify activity. "H 2 O" clumps are associated to one or more H 2 O maser from the BGPS GBT survey, Arcetri, and/or HOPS. Similarly, "CH 3 OH" clumps are associated to one or more CH 3 OH masers from Pestalozzi et al. (2005), Arecibo, and/or MMB. "UCHII" clumps are associated to the CORNISH UCHII regions. In section §5 we will use these categories to compare physical properties. Table 4 lists the detection statistics for sources matched to the BGPS in the overlap region and number of clumps matched to sources. We find that in the common overlap region 2223/4683 (47.5%) are starless candidates and 2460/4683 (52.5%) show some indicator of protostellar activity. This fraction of starless clump candidates is lower than the upper limit of 67% estimated from the ATLASGAL survey by Csengeri et al. (2014) which compare ATLASGAL clumps to WISE 22 µm and MSX 21 µm point source catalogs. The smaller fraction of starless clump candidates in our analysis is likely due to the inclusion of protostellar indicators with lower luminosity sensitivity such as Herschel 70µm sources in our analysis. The Herschel 70 µm sources accounts for the largest fraction of star formation indicators detected toward Protostellar clumps (2168/2460 or 88% of Protostellar clumps contain a 70 µm source flag). No strong variation in the Galactic longitude distribution is observed between Starless Clump Candidates and Protostellar Clumps ( Figure 5). Individual star formation indicator flags for each clump are listed in Table 5.
Except for 70 µm Unique, the BGPS indicator groups are not mutually exclusive. Indeed, 70/170 (41%) of UCHII clumps are also associated to a CH 3 OH maser and 103/170 (61%) are associated to a H 2 O maser. There are no BGPS clumps uniquely associated with H 2 O, CH 3 OH, or UCHII flags; those clumps are always associated with another protostellar indicator. This fact indicates that it is unlikely that a population of H 2 O unique clumps (meaning H 2 O masers are the only indicator of star formation activity) exist toward clumps that have been targeted by our GBT observations ( §2.2 and §3.3.1). Figure 6 shows the intersection between indicators with the fraction of the total calculated for each. For example, 70 clumps are associated to both CH 3 OH masers and UCHII regions: 0.412 (70/170) of those with UCHII regions are also associated to CH 3 OH masers, alternatively 0.235 (70/298) of those with CH 3 OH masers are also associated to UCHII regions.
Among these indicators the CH 3 OH and UCHII clumps are the most consistent indicators of high-mass protostellar activity. However, it is important to note that more extreme indicators of star formation activity do not necessarily represent a later stage in a linear clump evolutionary sequence. For one, not all clumps will have sufficient mass to form an OB star to drive an observable CH 3 OH maser or UCHII region. The sequence of star formation indicators outlined in Battersby et al. (2010) is consistent with the expected sequence of observational indicators for a high-mass protostellar core, but BGPS clumps are hierarchical constructs of density, likely composed of multiple cores (Merello et al. 2015). The surveys also do not have high enough resolution to identify single protostellar sources at large distances. For example, what appears to be a Hi-GAL compact source could be a protocluster at several kpc. In this case, the only certain true evolutionary transition is from the starless clump phase to the protostellar phase.

Comparison with Published Starless Clump Catalogs
We compare our protostellar indicator flags to the previously published starless clump catalog of Tackenberg et al. (2012) which identified 210 starless candidates with peak column density N(H 2 ) > 10 23 cm −2 from ATLASGAL 870 µm images in the 10 • < l < 20 • range. Tackenberg et al. (2012) use a combination of mid-IR colors from Spitzer imaging plus visual analysis of 24 µm point sources to identify starless candidates. Within this longitude range, there are 52 BGPS catalog sources with peak positions that lie within the BGPS angular resolution (33 ) of the Tackenberg et al. (2012) starless candidate position. We find that 30 out of those 52 clumps (58%) contain an indicator of star formation activity in our catalog. Nearly all of those associations result from identification of a 70 µm source in our catalog. This result highlights the necessity of far-IR imaging to identify deeply embedded and lower luminosity YSOs.
We also compare to the Traficante et al. (2015a) catalog of 1684 Hi-GAL clumps associated to IRDCs in the 10 • < < 55 • and |b| ≤ 1 • range. Clumps are extracted from the IRDC catalog of Peretto & Fuller (2009 Dunham et al. 2011b) and from morphological matching the BGPS clumps to the 13 CO spectra from the Galactic Ring Survey (GRS; see Ellsworth-Bowers et al. 2014). A Bayesian statistical methodology of formulating clump distance probability density functions (DPDFs) based on a range of prior distributions has been carried out and applied to the BGPS by Ellsworth-Bowers et al. (2013) and Ellsworth-Bowers et al. (2014). Further, to accurately propagate our uncertainty in heliocentric distance, we shall use Monte Carlo random sampling of the DPDFs when calculating the derived physical properties. The DPDF is calculated by multiplying the likelihood function, derived from the measured v LSR with a 7 km/s uncertainty to account for typical GMC dispersion and the Reid et al. (2014) rotation curve of the Galaxy, with prior distributions

Distance Resolution Broadcasting
The well-constrained BGPS distance sample contains 1414 clump DPDFs representing 43% of the BGPS velocity catalog in 10 • < < 65 • . Clumps show an easily distinguishable tendency to cluster in velocity coherent groups. Indeed, 75.4% of clumps in the overlap region have at least one adjacent neighbor (i.e., pixels touching) in the Bolocat label-masks. Clumps that are within spatial and kinematic proximity are more likely to be cospatial than chance alignments of noncospatial near/far kinematic distances. In order to increase the size of the Distance Sample, we apply a modified Friendsof-Friends algorithm (Huchra & Geller 1982) to associate or "broadcast" the known DPDFs based on whether clumps are nearby in ( , b, v LSR ) or position-position-velocity (PPV) space. Because this technique uses only proximity in PPVspace, the number of new associations does not have a strong selection effect from evolutionary stage or physical properties. To evaluate the optimal search angle-velocity search parameters we use the Distance Sample as a training set, and minimize the fraction of near/far misassociations.
We assign a clump to a "PPV-group" or "group" by searching a distance ρ around the clump in angular separation, θ s , and velocity, v s : Clumps within this neighborhood are assigned to a parent group and the group is iteratively expanded by searching around the neighbors. This algorithm effectively creates groups out of a minimum density in PPV-space. We calculate group assignments over a 30 × 30 grid of angular search and velocity search radii from θ s = 1.8 − 10.5 and v s = 1 − 6 km s −1 . We use the known maximum likelihood distances as a training set to evaluate the accuracy of the search parameters. We denote the conflict ratio, R conf , as the number of clumps in groups with disagreeing distance resolutions divided by the number of clumps in groups with at least two distance resolutions. With R conf as a figure of merit, we select the search parameters that maximizes the number of new distance resolutions for fixed R conf . It is important to note however that the DPDFs are not simply Boolean near/far kinematic distance resolutions. Thus, R conf can only be approximately as low as the mean weight between the near/far probabilities in the DPDFs. (Ellsworth-Bowers et al. 2013) find that 8% of the DPDFs disagree with the 13 CO Galactic Ring Survey kinematic distance ambiguity resolutions (KDARs) from Roman-Duval et al. (2009) using HI self-absorption (HISA). We fix the conflict ratio to be reasonably consistent with this figure at R conf = 0.12 and select the search parameters (θ s , v s ) = (3.9 , 4.0 km s −1 ) that maximize new associations. In total, we associate 446 new DPDFs over the full survey, and 226 in the survey overlap region. The KDA resolution broadcasting method thus increases the Distance Sample used in this study by 226 (16%) to 1640 clumps in the overlap region.
The DPDFs added to the Distance Sample are constructed from the weighted priors of the DPDFs within the group. The optimized search parameters are used to associate the neighboring clumps in the group to a common kinematic near/far distance. To create a new posterior DPDF for a clump in a group, we apply a weighted average of all priors from the clumps in the group with a DPDF to the the new group-member clump. The priors are given Gaussian weights by the PPVdistance ρ with FWHM values set to the optimized θ s and v s above. This weighting ensures that the closest clump priors in the group are most strongly weighted. To propagate our uncertainty in broadcasting the distance resolution of an individual clump, we apply a weighted error function (as a step-function analogue) to the new posterior DPDF to account for R conf = 0.12. Down-weighting the distance resolutions by the conflict fraction accounts for the mean uncertainty in the group associations.  Dunham et al. (2011b), 403 match to v2.0 sources and 100 have well-resolved KDAs in our sample. For sources flagged as either "near" or "tangent" in Dunham et al. (2011b), we find that 69/78 ≈ 88% agree. For sources flagged as "far" that had no incidence of an 8 µm absorption feature nor HISA, we find that only 11/20 ≈ 53% agree.

Distance Comparisons
The lack of an HISA feature is commonly used to associate GMCs to the far kinematic distance (Roman-Duval et al. 2010), where the superposition against the background HI distribution produces absorption at the velocity coincident with the highcolumn density present in the GMC. The clumps found in the BGPS are both on smaller spatial scales and at higher spatial densities than GMCs. The disagreement for far distance resolutions for clumps is because the HI distribution does not follow that of the clump dust continuum, and thus the lack of a HISA feature can not be used to accurately resolve the kinematic distance ambiguity for BGPS clumps. Wienen et al. (2015) use associations of 13 CO GRS isocontours in PPV-space and the associated HI absorption with those GRS clouds to resolve 44% of ATLASGAL clumps to the far kinematic distance. This is twice the 22% fraction found in Ellsworth-Bowers et al.

ANALYSIS OF PHYSICAL PROPERTIES WITH STAR
FORMATION INDICATORS In the following subsections, we shall compare and analyze observable and physical quantities for both the starless clump candidates and protostellar clumps. For the protostellar clumps, we sort clumps by their star formation indicators into subsamples described in §3.5. As noted in §3, this ordering is not an evolutionary sequence but is roughly ordered by the typical luminosity required to produce each indicator. In the following figures, the top two panels indicate the starless clump candidate and protostellar candidate subsamples while subsequent panels indicate various star formation indicator subsamples. For the physical properties and sampling techniques described below, Table 6 lists the derived physical properties and uncertainties for each clump in the Distance Sample, and Table 7 lists the statistical properties and uncertainties for subsamples of clumps sorted by star formation indicator.

Monte Carlo Sampling
For each subsample we compute Monte Carlo (MC) simulations to randomly sample all of the input parameters, both measured quantities with statistical uncertainties (e.g., S ν , T K ) and DPDFs. By calculating descriptive statistical quantities (e.g., µ 1/2 , ρ spear ) from suites of MC simulations, we can estimate the uncertainty and significance of comparisons between samples of clumps associated with different star formation indicators. The suites of MC simulations also more accurately propagate the uncertainties when compared to simply creating a single histogram of the maximum likelihood values, because some quantities such as heliocentric distance can have highly asymmetric and/or bimodal PDFs. Thus, while a single clump may have large fractional uncertainties in any given derived physical property, we can make robust and accurate calculations for large (∼ 10 2 − 10 3 ), statistically significant subsamples of clumps.

AGB Contamination Re-sampling
Contamination from AGB stars is a significant source of uncertainty when the only indicator of star formation activity within a clump is from the mid-infrared colors of a point source. In order to propagate this uncertainty due to AGB contamination, we re-sample flags discriminating between R08 YSOs and AGBs sources in the MC simulations. We use conservative estimates of the contamination fractions between the R08 YSO contaminants in the AGB classification and AGB contaminants in the YSO classification. As described in §3.2.2, we apply contamination fractions of 40% for AGB sources and 50% for YSO sources. All R08 sources within clumps are re-sampled this way. If a clump is associated to another star formation indicator besides R08, then it is still flagged as protostellar and the MC flag re-sampling will not affect the clump's overall designation. In the overlap region, only 280 clumps (6.0%) are uniquely associated with R08 sources, and 122 clumps (7.4%) in the Distance Sample. Resampling propagates our uncertainty in the R08 flags, but because of significant coincidence with other star formation indicators, it introduces only a minor effect in discriminating between starless clump candidates and protostellar clumps. In the following analysis, AGB contamination re-sampling is used in all MC calculations.

DPDF Sampling and Distance Biases
For the calculation of quantities which are distance dependent (e.g., R, M ), we draw the heliocentric distance from the source DPDF. Each MC simulation consists of 10 4 random draws. Since the output distributions of the physical properties may be skewed, we use robust statistical indicators such as the median (µ 1/2 ) and median absolute deviation (σ 1/2 ) 12 as conservative characterizations of the distribution of properties for a particular subsample.
Before we compare physical properties among clumps with different star formation indicators ( §5.2-5.4), we compare the distribution of DPDFs for each star formation indicator subsample defined in §3.5. Figure 8c shows the mean distribution or average DPDF for each indicator. The indicator categories show similar distributions of heliocentric distance. While the number of clumps associated with each indicator vary, the range and distribution of d are similar, with medians between 4.4 − 5.1 kpc. As a result, physical properties which depend linearly or quadratically on distance will not be strongly biased by the underlying distance distribution for each indicator subsample. Sources with associated CH 3 OH masers and/or UCHII regions show a relative deficit of distances from 2 − 4 kpc. This could be an effect from the method of resolving distances with the 8 µm absorption prior, where better spatially resolved clumps at the near distance are dominated by mid-IR emission from a high-mass YSO. These regions may not have the necessary contrast in the GLIMPSE 8 µm maps to result in a well constrained distance estimate, whereas further, and thus larger, clumps could have sufficient surrounding 8 µm absorption. In addition, the PPV method of associating nearby clumps with known distances for a fixed radius in angle and a fixed velocity better associates regions at larger distances where the angular separation between clumps is smaller.
When comparing distance independent quantities between starless candidate and protostellar categories, the total number of sources is nearly equal (2223 SCCs vs. 2460 protostar- containing clumps); however, when comparing distance dependent physical properties, it should be noted that in the Distance Sample there are 1.675 times as many protostellar-containing clumps than starless clump candidates (613 SCCs vs. 1027 protostellar clumps). The fraction of clumps in the Distance Sample that contain at least one star formation indicator is shown as a function of heliocentric distance in Figure 9. The median protostellar clump fraction in the Distance Sample is 0.60, and there is no significant trend in this ratio out to heliocentric distances of 12 kpc. Since only 58 clumps (3.5%) in this study have maximum likelihood distances greater than 12 kpc, this result indicates that differences in the properties of starless candidates compared to protostellar clumps are not driven by a bias in the protostellar fraction by distance.
We can only calculate distance-dependent quantities from the subset of clumps in each category of the Distance Sample. We run MC simulations for both the subsamples that have a well-constrained DPDF and those without to test how this selection criteria may bias our distributions. Figure 8a shows the observed total flux density S Total 1.1 for each indicator, with the total sample (light gray) and Distance Sample (dark gray). Figure 8b shows the cumulative results of 10 4 draws from the Monte Carlo simulation of the 1.1 mm flux density Gaussian-deviates and R08 AGB contamination re-sampling. The subsamples with DPDFs have greater median S 1.1 when compared to full samples. Similarly, figure 8d shows the distribution of the ratio of medians for each MC simulation. All indicator categories show a factor of 1.5 − 2.2 higher median flux for clumps with DPDFs; however, there is not a consistent trend in the ratio between star formation indicators. The ratio of median flux densities for protostar-containing clumps is only a factor of 1.1 − 1.2 times larger than the corresponding ratio for starless clump candidates. If we make the simple assumption that sources without well constrained DPDFs follow the average DPDF distributions for the Distance Sample, then this result implies that physical quantities that depend linearly

Flux Density, Flux Concentration, and Size
The observed 1.1 mm flux densities increase towards more extreme indicators from starless clump candidates to UCHII associated clumps (Fig. 8a). The median total 1.1 mm flux density of starless clump canddiates is 0.356 ± 0.005 Jy and is a factor of 15 lower in total flux density compared to the median for UCHII associated clumps. The starless candidate and UCHII associated clump distributions are distinct where 99% of starless candidates have lower total 1.1 mm flux density than the median total flux of UCHII associated clumps. Even the presence of a deeply embedded YSO observed as a compact Hi-GAL 70 µm source preferentially increases the median observed flux density by a factor of ∼ 1.8. Figure 9 shows these trends in an alternative way as a strongly increasing protostellar clump fraction with increasing 1.1 mm flux density. The ratio crosses 0.5 around 500 mJy and approaches 1.0 above 5 Jy. Without internal star formation activity, starless clumps would only be externally heated by the interstellar radiation field (ISRF; Evans et al. 2002, Jørgensen et al. 2006, Wilcock et al. 2012. As a result, clumps with weaker internal heating will have colder dust temperatures and lower observed 1.1 mm flux densities (see §5.3).
The deconvolved clump equivalent angular radius θ eq is calculated by subtracting the BGPS beam θ HPBW = 33 (Ω beam = 1234 sq. arcsec) in quadrature from the observed equivalent angular radius θ obs = Ω obs /π calculated directly from the sum of the 7.2 × 7.2 map pixels in the Bolocat label-mask frames where Ω obs is either the total solid angle defined by the Bolocat v2.0 seeded watershed algorithm or the FWHM solid angle defined by the contour at half-maximum flux density. Quantities that are derived from the Bolocat label-masks of the BGPS maps are sensitive to the spatial filtering imposed by the principal component analysis performed on the Bolocam time stream data. A clump size based on the FWHM definition is more robust to spatial filtering and is also more prevalent in the literature (e.g., Beuther et al. 2002;Shirley et al. 2003). In the following analyses, we shall use both total size definition, i.e., the full label masks, as well as the FWHM size definition to calculate derived properties. Figure 10a shows Ω Total for each indicator, with an increasing trend in the distributions with median 3.62 × 10 3 square arcsec (1.01 square arcmin) for the starless candidates to 1.81×10 4 square arcsec (5.03 square arcmin) for the UCHII associated clumps. There is also a systematically increasing trend in the concentration of the 1.1 mm flux density distributions when sorted by indicator. Figure 10b shows the ratio in deconvolved solid angles, or Total to FWHM "concentration", defined as Ω Total /Ω FWHM . Clumps that are near uniform brightness and with a peak brightness less than twice the 1σ (∼ 0.1 − 0.5 Jy) cut-off threshold of the clump boundaries should have concentrations near unity, while clumps with highly concentrated brightness distributions will have large concentration values ( 1). A strong increasing trend is seen from a median ratio Ω Total /Ω FWHM for the starless candidate phase at 1.55 to 16.38 for UCHII clumps. Clumps associated with H 2 O masers, CH 3 OH masers, and UCHII regions are concentrated at similar levels and also more concentrated than those with less extreme indicators.
The equivalent physical radius R eq is calculated from the observed Ω when projected to a heliocentric distance d drawn from the DPDFs R eq = 0.29 θ eq arcmin d kpc (4) Figure 11 shows the equivalent radius distributions for the Total (11a) and FWHM (11b) definitions. The median R Total eq increases from 0.842 ± 0.017 to 2.66 ± 0.08 pc from starless candidates to UCHII clumps; however, little separation among different indicators is observed using the FWHM definition with median R FWHM eq ∼ 0.6 pc. The R Total eq in this study are larger than those found by Schlingman et al. (2011) (median R eq = 0.75 ) for a subset of known 529 sources taken from the BGPS v1.0. This discrepancy can be attributed to the recovery of flux at lower spatial frequencies in the BGPS v2.0 maps resulting in larger clumps (Ginsburg et al. 2013). The different trends observed between the total source size and FWHM size are more striking and can be explained by offsetting factors. Clumps containing indicators of more luminous protostars (i.e., CH 3 OH masers, UCHII) have the highest flux concentrations which partially offset their larger total solid angle resulting in their FWHM sizes remain nearly the same on median as starless clump candidates and clumps with lower luminosity protostellar indicators. This result highlights the need to explore multiple size definitions in the analysis of source physical properties.

Gas Kinetic Temperature and Linewidth
NH 3 is an excellent intermediate to dense gas tracer with an effective excitation of ∼ 10 3 cm −3 for the lowest pure inversion transition (Shirley 2015). Fitting the intensities of the lowest energy inversion transitions provide estimates of the gas kinetic temperature (Ho & Townes 1983). The rotation temperature of the (1,1) and (2,2) p-NH 3 inversion transitions saturates for T K > 30 K, but the simultaneously observed (3,3) transitions retain sensitivity to T K < 100 K if the o-NH 3 (3,3) transition is not masing (Danby et al. 1988). Figure 12 shows the distributions of NH 3 -derived gas kinetic temperature with values ranging between T K ∼ 10 − 100 K. The starless candidate phase has a median T K = 13.96 ± 0.10 K with an increasing trend to T K = 27.2 ± 0.2 K for UCHII clumps. The clumps uniquely associated with compact Hi-GAL 70 µm sources, a sign of deeply embedded protostellar activity, show only a slight enhancement in the gas kinetic temperature. As a clump evolves from a quiescent phase to one with active star formation, radiative heating from star formation is expected to raise the gas kinetic temperature. With more extreme indicators of star formation activity suggesting higher luminosity YSOs, the CH 3 OH and UCHII region containing clumps are subject to stronger radiative heating from massive star formation. The starless clump candidate distribution does display a positive skew. This could be a result of undetected embedded sources or enhancement by the local radiation field from neighboring star forming regions. The lower S 1.1 in the candidate starless clump stage shown in Figure 8 can, in part, be attributed to this observed trend in T K . The flux density of a modified blackbody is 4.9 times higher at T K = 30 K than 10 K at 1.1 mm. We shall account for the temperature differences observed among the different star formation indicators in our calculations of the clump mass in §5.4.
The observed FWHM linewidths of NH 3 for clumps may be broadened due to a variety of factors including thermal broadening, bulk motion, optical depth, and turbulence. In the following analysis we use the NH 3 rather than HCO + linewidth because the NH 3 observations have three times better spectral resolution, 0.3 km s −1 , and NH 3 has lower optical depths, τ NH3(1,1) ∼ 1 − 5 compared to τ HCO + (3−2) ∼ 10 (see Shirley et al. 2013  served NH 3 linewidth, indicating a factor of two increase from 1.339 ± 0.014 km s −1 on median towards starless candidates to 2.736 ± 0.03 km s −1 on median towards UCHII associated clumps. For gas temperatures between 10 − 30 K the thermal linewidth is ∆v therm (NH 3 ) = 0.16 − 0.28 km s −1 . The observed NH 3 linewidths at > 1 km s −1 suggest supersonic turbulence exists on scales smaller than the 33 beam of the BGPS. The quiescent starless candidate phase and the deeply embedded candidate phase of uniquely identified Hi-GAL 70 µm compact sources have ∆v turb /∆v therm ∼ 10. The observed increase in linewidth with more extreme star formation indicators could be due to an increase in turbulent feedback in more luminous protostellar clumps. More luminous protostars will drive more powerful jets and outflows (see Bontemps et al. 1996) as well as provide more radiative feedback on the surrounding clump that can drive turbulence (see Krumholz & Thompson 2012). Alternatively, the larger linewidths towards clumps with more extreme indicators could be due to larger bulk-flow motions that are unresolved at the GBT angular resolution or a systematic tendency for higher mass stars to form in more turbulent regions. Without higher resolution observations of the clumps, it is difficult to determine which scenario drives the observed trend.
Prior observations of a significantly smaller subsample of BGPS clumps have indicated a breakdown in the linewidth-size relationship (Schlingman et al. 2011). We plot the linewidthsize relationship for SCCs and protostellar clumps observed in NH 3 , shown in Figure 13. The Spearman's rank correlation coefficients for SCCs and protostellar clumps are ρ sp = 0.24 and ρ sp = 0.50 respectively from the MC simulations. The starless clumps appear uncorrelated while protostellar clumps have a weak correlation that is only marginally better than the Schlingman et al. (2011) correlation (ρ sp = 0.40). SCCs, as traced by NH 3 emission, seems to have decoupled from the supersonic scaling relationship observed toward molecular clouds more so than protostellar clumps. Indeed the observed NH 3 model-fit velocity dispersion (expressed as FWHM linewidth) are below the extrapolated relationship observed by Larson (1981) for molecular clouds traced by CO 1 − 0. While these clumps themselves are still turbulent structures, their level of turbulence appears dissipated compared to the expected relationship in clouds when observed on clump spatial scales.

Mass Surface Density
The mean mass surface density is calculated using: where B 1.1 is the Planck function, κ 1.1 is the dust opacity, ζ is the dust-to-gas ratio, here assumed to be 1/100, and Ω is the source solid angle, and S 1.1 is the source integrated flux density. Both the Total and FWHM size definitions are used for Ω or S 1.1 . We assume the 1.1 mm emission is optically thin and use Gaussian-deviates to sample the statistical uncertainty. We also assume OH5 dust opacities calculated for coagulated grains with thin ice mantles (Ossenkopf & Henning 1994). For the Monte Carlo calculation, we draw Gaussian-deviates from the observational uncertainty on T K , and assume that the dust temperature is equal to the gas kinetic temperature. Because 47% (772/1640) of clumps with a distance lack T K measurements, for those clumps we draw T K from that indicator's observed distribution of T K (Fig. 12a). This takes into account the trend of increasing gas T K for clumps associated to more extreme star formation indicators. For clumps associated with multiple star formation indicators, the indicator with the hottest median T K is used. Figure 14 shows Σ Total H2 with an increasing trend towards higher Σ Total H2 for clumps associated to more extreme star formation indicators. The median Σ Total H2 increases from median values Σ Total H2 = 0.0145 ± 0.0002 g cm −2 and Σ FWHM H2 = 0.0165 ± 0.0003 g cm −2 for the starless candidate clumps to Σ Total H2 = 0.0207 ± 0.0005 g cm −2 and Σ FWHM H2 = 0.076 ± 0.002 g cm −2 for the UCHII associated clumps. A similar trend in increasing column density from quiescent (potentially starless) clumps to protstellar clumps has also been observed by other (sub)millimeter Galactic plane surveys (see Csengeri et al. 2014;Guzmán et al. 2015).
The factor of four increase in Σ FWHM

H2
is consistent with the greater flux densities and flux concentrations observed towards clumps with more extreme indicators. The median mass surface density for protostellar sources is close to the purported threshold for star formation of approximately 120 M pc −2 (or 0.025 g cm −2 ; Lada et al. 2010;Heiderman et al. 2010). Since half of these clumps presumably contain protostellar sources despite being below this threshold, this result highlights that BGPS clumps themselves contain significant substructure which are at higher mass surface densities.

Clump Mass
We calculate total clump masses with: when sorted by indicator from 10 4 MC simulations that account for R08 AGB catalog contamination and uncertainties in S 1.1 , T K (or randomly drawn from the appropriate T K distribution if no T K observations exists), and the DPDFs. The masses show a systematic increase from µ 1/2 (M H2 ) = 226 ± 11 M for SCCs and µ 1/2 (M H2 ) = 600 ± 20 M for protostellar clumps. Among the clumps associated to star formation indicators, the masses increase from µ 1/2 (M H2 ) = 390 ± 30 M for the 70 µm Unique category to µ 1/2 (M H2 ) = 2600 ± 200 M for the UCHII category. The difference in median masses between the SCC and Protostellar categories of ∆µ 1/2 (M H2 ) = 370 ± 20 M and the SCC and 70 µm Unique category of ∆µ 1/2 (M H2 ) = 170 ± 30 M . The median mass difference is driven by more SCCs than Protostellar clumps below 470 M and more Protostellar clumps than SCCs above that mass. We statistically estimate the highest mass SCC by integrating the upper-tail of M Total H2 PDF to where it equals one part in the sample size, 1 : 625, or the 99.84 th percentile. This sets the maximum observed mass of SCCs at < 1.4 × 10 4 M . The observed increase in mass towards clumps associated to more extreme star formation indicators could be due to several physical and systematic explanations. In the remainder of this section, we explore possible systematic effects to generate the observed mass difference.
For physical properties that depend on the dust temperature of the clumps, we assume an isothermal temperature with T dust = T K in the clumps. Nominally, these two temperatures are only expected to be well coupled by collisions at densities greater than ∼ 10 5 cm −3 (Goldsmith 2001;Young et al. 2004). Battersby et al. (2014) find no correlation between NH 3derived T K and Herschel-derived T dust for a quiescent clump; however, agreement within 20% is observed when the Herschel data is not background subtracted. Battersby et al. (2014) suggest that the gas and dust are weakly coupled at densities 10 4 − 10 5 cm −3 . Similarly, the T dust derived from far-infrared SEDs in Traficante et al. (2015a) show a 20% agreement to T K for both starless candidate and protostellar clumps. A 20% decrease in T dust (from T K ∼ 12) will increase our mass and column density measurements for the starless candidate group by 30%. This effect alone is not enough to explain the observed mass differences. There is no reason to decrease only the dust temperatures of starless clump candidates. Any simultaneous decrease in T dust for protostellar clumps would mitigate the decrease of the median mass differences.
In reality, there are dust temperature gradients within the BGPS beam due to heating from the interstellar radiation field and due to embedded protostars in protostar-containing clumps. The general effect of line-of-sight dust temperature gradients on greybody SED fitting of dust continuum emission is to overestimate the dust temperature and underestimate the mass more severely for starless clumps than protostellar clumps (Malinen et al. 2011). The corresponding effect on the line-of-sight average NH 3 T K cannot explain the observed mass difference unless the temperature gradients cause an underestimate of the average line-of-sight T K for the starless clump candidates from the observed median value of 12 K down to 8.8 K, a value that is lower than all NH 3 -based T K measurements in this paper and 95% of SED-based T dust measurements in Traficante et al. (2015a). Analysis of the dust temperature profiles from grids of radiative transfer models of low-mass starless cores (Shirley et al. 2005;Launhardt et al. 2013) find that the mass-weighted average T dust never drops below 10 K for a grid of Bonner-Ebert spheres (Ebert 1955;Bonnor 1956) with central densities spanning 10 4 to 3 × 10 6 cm −3 and strength of the interstellar radiation field equal to half the standard value (Habing 1968). These low-mass starless core models are not the exact analogues of massive starless clumps which are possibly fragmented and clumpy, have lower observed average density (Figure 14), and are likely subject to stronger heating from the interstellar radiation field than half the Habing value, effects which will increase the average T dust compared to low-mass starless cores. These results indicate it is unlikely that dust temperature gradients result in an average T dust as low as 8.8 K. Radiative transfer modeling of the entire BGPS clump population is beyond the scope of the current paper and requires higher spatial resolution continuum and NH 3 images. In the absence of reliable dust continuum temperatures for the full sample, a single component T dust = T K assumption is the next best alternative.
The other assumption made in calculation of mass surface density and mass is that the dust opacities are well described by OH5 dust. They are a popular set of opacities because grains are expected to coagulate and accrete ice mantle in dense environments and there are observational constraints that indicate OH5 opacities are a reasonable match, albeit toward nearby, low-mass cores (see Shirley et al. 2011). The OH5 model assumes that the grains have coagulated for a period of 10 6 years at a density of 10 5 cm −3 . Both this time period and the density used in the coagulation model are likely larger than the typical values observed toward BGPS clumps (see Schlingman et al. 2011;Dunham et al. 2011a). In order to explain the ob-   Krumholz & McKee (2008). The dash-dotted line shows the 0.025 g cm −3 threshold for star formation observed in local clouds (Lada et al. 2010;Heiderman et al. 2010). served ∆M Total H2 the opacity ratio towards protostellar clumps versus starless clumps would have to be 1.7 − 2.6 (for 70µm unique clumps or the full protostellar clump sample). While it is certainly likely that there are opacity variations among clumps that deviate from the assumed OH5 opacities, the variation would have to be systematic between starless candidate and protostellar clumps. There is currently no evidence for or against such a systematic variation in opacities.
The observed mass difference cannot be due to selection effects on the flux density S 1.1 of the Distance Sample alone (i.e., the subsample of clumps that have DPDFs compared to those without) because SCCs and protostellar clumps have similar ratios in median S 1.1 (see §5, Fig. 8d). At most this effect could account for ∼ 10%, lowering the mass difference to ∆µ 1/2 (M H2 ) ≈ 330 M . In addition, as a flux limited survey, the BGPS is effected by Malmquist bias and mass incompleteness. Ellsworth-Bowers et al. (2015)  = 230 M , is at the 80% completeness level, and the 90% completeness level is achieved at the 64 th mass percentile. The incompleteness in the SCC category should not, qualitatively, affect the observed median mass difference. Because R proto is smaller at low flux densities and masses, the population of undetected clumps with M < 400 would disproportionately add to the SCC category and simply enhance the median mass difference further.
One additional possibility is the spatial filtering in the BGPS survey has systematically affected the flux densities of sources. Ginsburg et al. (2013) estimate the effects of spatial filtering by approximating the source brightness distribution as a powerlaw with varying index and finds that the flux density for sources with angular radii less than 40 is essentially recovered (see Ginsburg et al. (2013) Figure 8). Depending on the exact value of the power-law index, the spatial filtering recovers > 75 − 80% of the flux density for sources with angular radii < 60 . This accounts for 77% (3606/4683) of the BGPS sources. It is evident from Fig. 10a that the sources more strongly affected by spatially filtering are protostellar stages with more luminous evolutionary indicators which tend to have larger angular radii. Thus, accounting for spatial filtering should systematically increase the median mass difference between starless clump candidates and protostellar clumps. A detailed analysis of the effects of spatial filtering on a source by source basis is beyond the scope of the present work.
In summary, the observed mass difference cannot be explained by biases between the Distance sample and the full sample, incompleteness, spatial filtering, or our assumptions in the mass calculation unless there is a factor 1.7−2.6 systematic variation in dust opacities between starless clump candidates and protostellar clumps. Interestingly, the BGPS sample is the second statistically significant sample to see this systematic offset in mass. The Herschel survey of (Traficante et al. 2015a) toward IRDCs also find a modest offset of ∼ 80 M between the average mass of starless clumps and protostellar clumps. While the mass difference is smaller by a factor of two in their study than ours, the mass difference in their survey is significant. Alternatively, Hoq et al. (2013) using data from the first year data release of the MALT90 molecular line mapping survey of several hundred ATLASGAL clumps in the fourth and first quadrants do not observe an increasing trend in total mass from an "8/24 µm quiescent" phase to a "PDR" phase. It should be noted that these two surveys are drawn from subsets of clumps in the Galactic plane; therefore interpretation of these results depend on how the properties of the subsets compare to the more complete clump populations observed in (sub)millimeter continuum Galactic plane surveys. In §6.3 we explore possible physical explanations for the variation in the median mass between the starless clump candidates and protostellar candidates under the assumption that the observed mass difference is real.

Virial Parameter
The virial parameter, α = M vir /M , expresses the relative importance of gravitational and kinetic energy, and can be used to assess whether clumps are gravitationally bound and/or where a 1 = 1−p/3 1−2p/5 is the correction factor for a power-law density distribution with index p and a 2 ∼ 1 for sources with aspect ratios less than 2 (Bertoldi & McKee 1992). For a clump with negligible magnetic fields, α = 1 is gravitational virial equilibrium and α ≈ 2 is marginally gravitationally bound (see Kauffmann et al. 2013). The BGPS does not have sufficient resolution to constrain p, so lacking appropriate data, we draw p = 1.8 ± 0.4 Gaussian-deviates in the MC simulations, as constrained from dust continuum modelling of high-mass star forming regions . Here we use the total mass and radius, thus setting the outer boundary to calculate α at the extent of the 1.1 mm emission. Figure 16 shows the virial parameter calculated from MC simulations for starless candidates and protostellar clumps with median values α = 0.73 ± 0.06 and α = 0.68 ± 0.03 respectively, suggesting that most clumps are in approximate virial equilibrium with more than half of clumps in the Distance Sample with NH 3 observations showing (sub)-virialized motions (α 1). Further, 76% of the starless candidates and 86% of the protostellar clumps have α < 2 implying that most BGPS clumps are gravitationally bound. Although 24% of the starless candidates are unbound by the criteria α > 2, it is probable that these clumps host gravitationally bound higher density substructures. These results suggest that the majority of BGPS clumps in the Distance Sample are self-gravitating or collapsing.
Observations of GMCs indicate virial parameters that are typically 2 (Larson 1981;Solomon et al. 1987;Scoville et al. 1987). For instance, analysis of GMCs in the Galactic Ring Survey find a mean value of α = 1.9 ). For smaller clump scales there exist a range of virial parameter measurements that report samples with α 1 or α 1 (see Kauffmann et al. 2013 for a summary of literature values).
We find our measurements of α 1 consistent with other Galactic plane surveys towards clumps from both BGPS and ATLASGAL with NH 3 -determined linewidths. The most direct comparison are previous studies using GBT observations of NH 3 that have been carried out towards the BGPS sources (Dunham et al. 2010(Dunham et al. , 2011b. Dunham et al. (2011b) measured the virial parameter towards 456 BGPS clumps with a median α = 0.74. Wienen et al. (2012) also measure the virial parameter toward ATLASGAL clumps using NH 3 , and find a mean α 1 for clumps.
One important aspect that is often ignored in the standard virial parameter analysis is the importance of magnetic fields in virial balance. Magnetic fields in dense gas are notoriously difficult to measure, and there have been no systematic observations toward BGPS clumps. We can estimate the B-field strength required for virial balance from the expressions derived in Kauffmann et al. (2013)  The required B-field strength scales as B ∼ 1/α, implying that very sub-virial clumps with α 0.1 required B 600 µG. In these more extreme cases, B-field support may prove difficult, but for the typical (median property) starless clump candidates, the required B-field strength for support is within a reasonable range of observed values (for a summary of B-field measurements, see Crutcher 2012).
6. DISCUSSION For many of the observed properties and calculated physical properties analyzed in §5, there is a systematic trend in the property with the ordering of the star formation indicators from 70 µm unique flags, mid-IR star formation flags, H 2 O maser flags, CH 3 OH flags, to UCHII flags. Protostellar clumps that contain the most extreme indicators of luminous protostars (CH 3 OH maser and UCHII containing clumps) tend to have the highest flux densities, are the most centrally condensed, have the highest temperatures, are more turbulent and more massive on median than other clumps. A partial evolutionary analysis of the ATLASGAL clumps associated to RMS sources, CH 3 OH maser, and UCHII regions, has been performed (Urquhart et al. , 2013a(Urquhart et al. ,b, 2014. The BGPS results are consistent with the observed trends (M , T K ) for the properties of ATLASGAL clumps that contain methanol masers and HII regions (see Table 4 of Urquhart et al. 2014). In contrast, the observable and physical properties indicate that starless clump candidates have lower flux densities, are less centrally concentrated, have smaller sizes, are colder and less turbulent, and are less massive on median when compared to clumps with indications of protostellar activity. In this section, we explore the potential for forming massive stars, the timescales for evolution of this newly discovered population of starless clumps, and explore possible explanations for the observed median mass difference between starless and protostellar clumps.

Massive Star-forming Potential of Clumps
With a robust sample of starless clump candidates, it is important to assess their potential to form high-mass stars. We can estimate the mass of the most massive star formed by a starless clump candidate, M max , using the stellar Initial Mass Function (IMF) from Kroupa (2001) 13 . The total stellar cluster mass is equal to the mass of the progenitor clump times a star formation efficiency factor SF . We can derive a relationship between M max and the mass of the clump from   Figure 17 plots the number of clumps in the Milky Way capable of forming stars with masses ≤ M max . Note that the number of clumps is very sensitive to assumed star formation efficiency, especially for very massive stars.
An alternative way to assess the massive star-forming potential of BGPS clumps is to compare M H2 and R eq ( Figure  18). Values are represented as the median value drawn from the MC simulations per clump. A robust increasing correlation is observed with ρ sp = 0.942. To show the differences in clumps by star formation activity we plot SCCs in yellow and clumps associated to UCHII or CH 3 OH masers as unambiguous detections of high-mass star formation activity in red. The distributions of points follow the general trends in discussed in §5: SCCs are generally lower mass, smaller radius, and lower mass surface density. Several thresholds for star formation 14 The correction factor is calculated numerically by integrating the H 2 density distribution over the volume of the wedge observed by the BGPS divided by the integral of the H 2 density distribution over the total volume of the Galaxy. The Wolfire et al. (2003) density distribution is given by ρ(Rg, zg) ∝ exp(−4 ln 2(Rg − 0.571R ) 2 /(0.52R ) 2 ) exp(−|zg|/0.059) for R < 0.82R and ρ(Rg, zg) ∝ exp(−Rg/(0.34R )) exp(−|zg|/0.059) for R ≥ 0.82R where Rg is the Galactocentric radius in kpc, zg is the height above the Galactic mid-plane in kpc, and R = 8.5 kpc is the distance to the Galactic center. Note that we account for the vertical offset of 25 pc of the Sun above the Galactic plane in the calculation of the BGPS wedge integrals (see Appendix C   determined from local clouds are over-plotted. Fitting star formation rates to local molecular clouds, Lada et al. (2010) and Heiderman et al. (2010) find relationships for "effecient" star formation of ≈ 116 M pc −2 and 129± M pc −2 , respectively, where above these thresholds the star formation rate is linearly proportional to mass surface density. The thresholds from Lada et al. (2010) and Heiderman et al. (2010) approximately bisect the samples, with 50.6% and 42.5% of the full samples, 46.6% and 39.9% of the UCHII and CH 3 OH maser sample, and 41.4% and 33.6% of the SCC sample each exceeding the respective limits. Kauffmann & Pillai (2010) derive the more restrictive criteria of M ≤ 580M (R eq /pc) 1.33 from observations of nearby molecular clouds. Note that Kauffmann & Pillai (2010) scale their OH5 dust opacities down by a factor of 1.5, so for consistency to values in this work we use the leading factor of 580 M rather than the original 870 M . In the BGPS Distance Sample, 24.8% of the full sample, 54.9% of the UCHII and CH 3 OH maser sample, and 10.8% of the SCC sample each exceed this criteria. Because nearly half of BGPS clumps with an unambiguous detection of a high-mass YSO do not meet this criteria, this suggests that the BGPS cloud structures do not strictly follow the criteria derived in Kauffmann & Pillai (2010). While the majority of SCCs lie below this threshold, any systematic growth would increase clump total mass and radius, placing them in phase-space for a higher probability of high-mass star formation (see §6.3.2). In a similar blind sample from the ATLASGAL survey, Wienen et al. (2015) find that 92% of ATLASGAL clumps meet this criteria and 100% above the Lada et al. (2010) and Heiderman et al. (2010) lines. Notably, only 1 clump in the BGPS from 10 • < < 65 • (W49, G043.167+00.011) meets the criteria for the progenitors of Young Massive Clusters (YMCs; Longmore et al. 2014) with clump masses 10 5 M , although W49 is not a SCC because it is actively forming stars.
6.2. Clump Timescales 6.2.1. Average Free-fall Time For the SCC sample we estimate the average free-fall time, t ff , calculated with where n c is the clump central density. We estimate the clump central density n c from the peak mass surface density and FWHM radius, assuming a spherically symmetric Gaussian density distribution (cf. Appendix C of Pattle et al. 2015). Integrating over the FWHM cylindrical aperture yields the expression where θ is the deconvolved angular radius for the Total or FWHM definition (see §5.2) and µ is the mean molecular weight. While a density distribution described by a singly peaked Gaussian is a crude assumption, it is more representative of higher density sub-structures within starless clumps than just the average density with the FWHM. Figure 19 shows the average free fall time plotted versus the mass of the clumps of the MC simulation. There is no significant correlation between the average free-fall time with the mass of the clumps. The median is t ff,c = 0.84 ± 0.02 Myr for the starless clump candidates.
Since the BGPS clump FWHM sizes are typically a factor of two larger when compared to the Jeans length of 0.36 pc (T K /15 K) 1/2 (10 3 cm −3 /n c ) 1/2 , these starless clumps should fragment into smaller structures. Higher resolution observations of a subset of BGPS clumps show this fragmentation into objects with higher densities (Merello et al. 2015). As a result, the median t ff is likely an upper limit to the free-fall time for a typical starless clump candidate to form a single detectable star.

Clump Lifetimes from Number Statistics
If the Galactic population of clumps is in steady state where the clump formation rate is equal to the destruction rate then the number fraction of clumps in a given state is proportional to the lifetime of that state. We find nearly equal numbers of starless clump candidates (47.5%) compared to protostellar clumps. This result implies that the typical lifetime of a starless clump candidate is nearly equal to the typical lifetime of a protostellar clump. The starless clump lifetime is defined as the time it takes for a clump detectable by the BGPS to form a single star that is detectable in current infrared surveys of the Galactic plane. The protostellar clump lifetime is determined by the time it takes a cluster of stars to form (not a single star) and dissipate the gas and dust to the point where the clump is no longer detectable in the BGPS. This latter timescale is difficult to determine, but we can use the cluster formation timescale as a rough estimate. This cluster formation timescale has a large range and is typically longer than 1 Myr and up to a few Myrs (Rebull et al. 2007, Chomiuk & Povich 2011, Morales et al. 2013; see review by Longmore et al. 2014). The total population statistics imply a long average timescale for starless clump candidates that is on the order of a few times the average free-fall time of the clumps and that is consistent with the effects of turbulence and or magnetic fields lengthening the starless evolution timescale.
A single timescale cannot properly describe the starless clump evolution timescale because clumps with different masses should evolve out of the starless phase at different rates. While a median mass starless clump candidate (230 M ) may have a lifetime over 1 Myr, it seems highly unlikely that 10 4 M starless clumps exist in a starless phase for that long. The lack of very massive (M > 1.4 × 10 4 M ) starless clump candidates argues for a shorter timescale for these objects to form a single detectable star and move to the Protostellar clump category. Current single timescale estimates (or upper limits) in the literature for massive starless clumps (Csengeri et al. 2014;Traficante et al. 2015a) are problematic because they do not account for this variation with mass and because they also attempt to compare to a single timescale for the lifetime of the embedded massive star formation phase (i.e., Davies et al. 2011;Duarte-Cabral et al. 2013). For the remaining discussion, we analyze how starless clump lifetimes vary with the mass of the clumps.
We can estimate how the clump lifetime for massive clumps vary with mass by using a method pioneered in Battersby (2013, Ph.D. Thesis) by pinning the relative lifetimes of clumps to an estimated absolute lifetime of the Class II CH 3 OH maser τ CH3OH = 2.5 − 4.5 × 10 4 yr (van der Walt 2005). This timescale is a statistical estimate based on correcting the number of the observed CH 3 OH masers in Pestalozzi et al. (2005) for completeness to predict the total Galactic count of highmass stars with mass > 20 M . van der Walt (2005) uses a SFR = 4 M yr −1 based on the SFR = 3−6 M yr −1 range in Boissier & Prantzos (1999). Chomiuk & Povich (2011) re-analyze several Galactic SFR measurements (see Table 1 therein) to derive SFR = 1.9 ± 0.4 M yr −1 . The Chomiuk & Povich (2011) study is unique in that they apply the same IMF (Kroupa 2001) and SFR law to all measurements and use the Starburst99 code to self-consistently compare studies with substantially different methodologies (including free-free radio continuum observations as well as YSO counts). Scaling the CH 3 OH maser lifetime inversely by their revised SFR yields τ CH3OH = 6.1 − 9.3 × 10 4 yr.
The clump lifetime is calculated using the relative number fraction of sources compared to the CH 3 OH maser sample at a given mass via where τ CH3OH is the completeness corrected and the SFR corrected Class II CH 3 OH maser lifetime. The ratio of N SCC /N CH3OH is corrected for the fraction of SCCs and clumps with CH 3 OH that are in the Distance sample compared to the total sample. The fundamental assumption of this method is that all clumps capable of forming a 20 M star (M clump > 1000 M ) will at some time form high-mass YSOs with detectable CH 3 OH masers. A population of equal mass clumps that are not capable of producing observable CH 3 OH masers would decrease the observed fraction and lead us to overestimate the inferred lifetime. This method also assumes the lifetime of a single CH 3 OH maser (τ CH3OH ) for all clumps with CH 3 OH maser associations. For clumps associated to multiple CH 3 OH maser sites (i.e., YSOs) an age spread within the clump would increase the observed duration of the clump CH 3 OH maser phase and lead us to underestimate the phase lifetime when the shorter, single YSO τ CH3OH is applied; however, for clumps associated to MMB sources, which have accurate interferometric positions, only 25/272 (9.2%) are host to more than one CH 3 OH maser, so this does not likely introduce a strong bias. To calculate the lifetimes as a function of mass we extend this method to the relative fractions of MC samples in narrow mass bins (0.0125 dex logarithmic bins for 10 4 MC mass samples per clump). In effect this is the ratio of two mass PDFs each scaled by the sample sizes and then multiplied by τ CH3OH . Figure 20 shows the lifetime of SCCs as a function of mass. The starless candidate clump lifetime approximately follows τ SCC ∼ 0.37 ± 0.08 Myr 10 3 M /M for M > 10 3 M . The uncertainty corresponds the systematic uncertainty in τ CH3OH . It is not advisable to continue our timescale analysis to < 1000 M as the assumption that all clumps will form a CH 3 OH maser is less likely to be true; however, extrapolation of the trend to lower masses indicates that the starless clump lifetime becomes longer than the average clump free-fall time at < 450 M . This is approximately twice the median SCC mass. Thus, most starless clump candidates have lifetimes longer than their average free-fall times of 0.8 Myr. In turn, very massive starless clump candidates with M H2 > 10 4 M have very short lifetimes τ 0.03 Myr consistent with lack of such massive objects detected in the BGPS survey . Our lifetimes estimates are 5 times longer and 3.5 times longer than the lifetimes estimates of Tackenberg et al.  degree survey area to the entire Galaxy.
6.3. Origin of Starless vs. Protostellar Clump Mass Difference In §5.4.2 a systematic difference between the median mass of protostellar clumps and starless clump candidates of 170 − 370 M was discussed. Fundamentally, the median mass difference is due to more starless clump candidates with masses below M 470 M and fewer higher mass starless clump candidates than protostellar clumps (see Figure 21). We explore three possible physical explanations for this mass difference: (i) not all starless clump candidates will evolve into detectable protostellar clumps, (ii) clumps grow in mass by accreting surrounding material from their parent GMC, and (iii) clumps evolve at different rates from the starless to protostellar phase in exact proportion to their number statistics.

Infertile Starless Clumps
Perhaps the simplest physical explanation for the systematic mass difference is that there is a large population of clumps among the SCCs that, while starless, will not evolve to form YSOs detectable by our star formation indicators. In order for this population of objects to remain starless clump candidates, they would either remain "infertile" and never form stars or would only be able to form stars with bolometric luminosities that are less than the 30 − 140 L sensitivity (corresponding to 1 − 2 M protostars accreting at typical rates) of current far-infrared Galactic plane surveys ( §3). Such a population of low-density and low-mass clumps would bias the SCC distribution to lower masses. For a Kroupa IMF, a 230 M median mass SCC is capable of forming stars up to 6 M , which is detectable given the sensitivity limits of current far-infrared surveys (see Figure 4); however, for SCCs at the 50% BGPS mass completeness limit of 70 M , the maximum mass star formed is only 2.5 M which is closer to the sensitivity limit of current mid-and far-infrared surveys.
One test of this hypothesis is to check for the presence of dense gas. The HCO + 3 − 2 transition has an effective excitation density n eff > 10 4 cm −3 for the typical properties of starless clump candidates (Shirley 2015). HCO + 3 − 2 from Shirley et al. (2013) was detected towards 44% of SCCs and 78% of protostellar clumps in the Distance Sample. This result indicates that there is either a population of lower density SCCs or that their filling fraction of dense gas is less than toward protostellar clumps. Since the 3 − 2 transition has an upper energy level that is 25.7 K above ground, the effective excitation density (n eff ) is a very sensitive function of T K at the median temperature of 12.8 K measured toward SCCs. The sensitivity of the 3 − 2 transition to T K means that the detection is biased against the SCCs that are colder on average than protostellar clumps. So this test is biased meaning that since SCCs are colder on average than protostellar clumps, HCO + 3 − 2 requires higher densities in SCCs (by a factor of ≈ 2) to be excited than toward protostellar clumps.
Another test of this hypothesis is to carefully control systematic effects from sampling clumps with different physical properties by performing a variety of astrophysical cuts on the sample and calculate the distribution of ∆µ 1/2 (M H2 ) from the MC simulations. These cuts are designed to remove low-mass, low-density, or unbound objects. Figure 22 shows the cumulative distribution functions (CDFs) for ∆µ 1/2 (M H2 ) subsamples selected on mass, NH 3 -derived gas kinetic temperature, heliocentric distance, virial parameter, mass surface density, central density, and molecular detections. The observed mass differences between categories, ∆µ 1/2 (M H2 ) are robust to these criteria, with the values calculated with respect to the Hi-GAL 70 µm Unique category at ∆µ 1/2 (M H2 ) ∼ 100 − 200 M and the values calculated with respect to the Protostellar category ranging between ∆µ 1/2 (M H2 ) ∼ 200 − 500 M . The robustness of the median mass difference to astrophysical cuts does not support the explanation that a population of low star formation potential SCCs solely drive the increasing trend in clump total mass.

Possible Clump Mass Growth
BGPS clumps commonly reside in a complex structural hierarchy of clouds and filaments. When sorted by protostellar indicator, the increasing trend in both Σ and M tot but weak increasing trend in R eq are consistent with the physical process of clumps, as focusing points of local gravitational potential minima, accreting mass from the surrounding cloud in a "conveyor belt" fashion while maintaining the same approximate volume. Clump mass growth via accretion from the surrounding molecular cloud is a possible mechanism for the difference between SCCs and protostellar clump masses. An upper-limit to the clump accretion rate can be set by the free-fall time of the highest density sub-component of the clump. An estimate of t ff from the MC simulation shown in Figure 19 (2). The values ∆µ 1/2 (M H 2 ) ∼ 100 M agree across a wide range of physical cuts designed to remove low mass, low density, or unbound objects. explanation for the mass difference. Note that all molecular clouds cannot be collapsing on free-fall timescales and still be consistent with the observed Milky Way SFR (Zuckerman & Evans 1974). Other mechanisms such as magnetic and turbulent support would increase the clump collapse time and decrease the predicted growth rate. As a simple estimate, we calculate the total mass available to a clump embedded in a GMC with typical properties. We assume a GMC with a spherical inflow velocity of ∆v ∼ 1 km s −1 onto a clump, feeding from a zone of ∆R feed ∼ 1.0 pc over t ∼ 1 Myr. For median sized clump with R eq ∼ 1 pc and feeding zone density n ∼ 150 cm −3 (N.B. median n = 230 cm −3 in GMCs in the Galactic Ring Survey, Roman-Duval et al. 2009), the homogeneous, spherical-shell reservoir mass would be M res ∼ 250 M . Higher inflow velocity or density would both increase the available mass. These values show that reasonable physical conditions within a GMC yield estimates that are consistent with the observed median mass difference between the SCCs and protostellar clumps, suggesting that cloud-to-clump accretion as a plausible physical mechanism.
There are theoretical models of such large scale clump ac-cretion. For instance, Vázquez-Semadeni et al. (2009) invoke a model of hierarchical gravitational fragmentation in which multi-scale collapse occurs. In this picture, the mass accretion rate onto the clumps is determined by the clump tidal radius and initial gravitational well hierarchy (see also Smith et al. 2009). In another approach, Murray & Chang (2012) have explored a model of Bondi-Hoyle accretion of clumps within GMCs and predict that the clump mass accretion rate is slightly super-linear with mass of the clump (Ṁ ∼ M 5/4 clump ). Maschberger et al. (2014) instead use a model of stochastic growth with a sub-linear mass growth rate (Ṁ ∼ M 2/3 sink ). The exact nature of clump accretion and how it depends on the mass of the clump is still under debate; nevertheless, there is theoretical support for cloud-to-clump accretion.
There is also direct observational evidence for large scale flows onto clumps and filaments. One of the most striking examples is a flow estimated at 2500 M Myr −1 traced by HCO + 1 − 0 blue asymmetric line profiles along a filamentary complex into the central clump of the source SDC335.579 − 0.272 (SDC335; Peretto et al. 2013). The total mass of this region is 5500 M which would place this among the most massive clumps in our survey. Another example of inflow onto clumps is seen in the converging flows toward the DR21 filament complex with inflow rates of 1000 M Myr −1 onto 4900 and 3300 M clumps from infall profiles of the 1 − 0 transitions of HCO + and CO . Analysis of velocity centroids can also reveal potential mass inflow onto clumps along filaments (see three examples from the IRDC survey in Tackenberg et al. 2014) although there can be a degeneracy between inflow and outflow in the interpretation (Henshaw et al. 2014). These regions are already forming protostars and display hub-filament geometries where filaments of dense gas appear to converge on the central hub, so they are not exact analogues to the massive starless clumps in the BGPS survey; nevertheless, they show that clumps are deep gravitational potentials that draw surrounding material in.
The examples of large scale inflow also extend to nearby lower mass regions. For example, inflow signatures are seen in line profiles of HNC 1 − 0 toward the Serpens South region indicating mass flow rates of 130 M Myr −1 onto the filament complex and 28 M Myr −1 onto the central cluster (Kirk et al. 2013; for a higher resolution study of the same region, see Fernández-López et al. 2014). The central filament (L1495, B213) of the low-mass star-forming Taurus complex also shows direct evidence of a large scale flow in CO 1 − 0 of 27 − 50 M pc −1 Myr −1 perpendicular onto the filament complex (Palmeirim et al. 2013). While these last two examples occur in lower mass regions than the typical BGPS starless clumps, they show that flow rates of 10 − 100 M Myr −1 are possible in low-mass regions. It therefore seems plausible that higher mass inflow rates of a few 100 M Myr −1 should be possible in the more massive starless clump regions see in the BGPS.

Different Starless Clump vs. Protostellar Clump Lifetimes
A third explanation for the observed mass difference is that the lifetime of the SCC phase decreases with increasing total clump mass in exact proportion to the fraction of SCCs at each mass. The shorter lifetime of massive clumps explains the lack of very massive (M > 10 4 M ) starless clumps candidates. We estimate the lifetime of this phase is less than 0.03 Myr based on the CH 3 OH maser technique. Thus protostellar clumps with > 10 4 M require an alternative formation mechanism: the average mass-growth from cloud-to-clump accretion is not sufficient to explain the lack of these massive SCCs. Using even a high mass infall rate of ∼ 10 3 M Myr −1 , the most massive clumps could only accrete a meager few ∼ 10 1 M over their few ∼ 10 4 yr starless-phase lifetimes.
The lifetime ratio of SCCs to Protostellar clumps would have to decrease with increasing mass in proportion to the ratio N SCC /N proto if lifetimes arguments are the sole explanation for the observed median mass difference. The ratio of N SCC /N proto scales as M −0.4 over the range 100 ≤ M ≤ 1000 M . The negative exponent implies that the lifetime scaling for SCCs is a steeper function of mass than the Protostellar clump lifetime. The value of the exponent is insensitive to the lower bound of the mass between 100 and 300 M . There is no known theoretical prediction of such a lifetime scaling between high-mass SCCs and high-mass Protostellar clumps.
With the current data, it is difficult to disentangle the relative importance of differing lifetimes for SCCs and Protostellar clumps versus the importance of clump mass growth in explaining the median mass difference between SCCs and Protostellar clumps. It is likely that both processes occur.

SUMMARY
We present a sample of 4683 molecular cloud clumps from the BGPS sorted by different observational indicators of star formation activity. We use a variety of Galactic plane surveys in a common overlap region between 10 • < < 65 • that include: 70 µm compact sources from Herschel Hi-GAL, mid-IR color-selected YSOs, H 2 O and CH 3 OH masers, and UCHII regions. We use Monte Carlo random sampling to calculate the clump physical properties using a subsample of 1640 wellconstrained clump DPDFs. We also present a catalog of 1663 clump NH 3 T K measurements and 22 GHz H 2 O maser observations.
We find the following conclusions: 1. We identify a subsample of 2223 dense clumps with no indicators of star formation activity, representing the largest and most robust sample of starless clump candidates from a blind survey to date.
2. We measure numerous increasing trends in median physical properties from starless candidates to the most extreme indicator of star formation, UCHII regions. These include: S Total 5. Virial parameters derived from NH 3 show sub-virial clumps with median α = 0.7. More than 75% of BGPS clumps are gravitationally bound (α < 2). A medianproperty SCC would be in virial equilibrium if a modest ∼ 50 µG magnetic field is present.
6. We find a median mass difference between the SCC and prostellar categories of ∆M Total H2 = 170 − 370 M . It is unlikely that this mass difference can be solely explained by a population of SCCs which are incapable of forming detectable protostars in current Galactic plane surveys. If the observed mass difference is not due to a systematic increase in dust opacity of 1.7 − 2.5 from SCCs to protostellar clumps, then there are two likely explanations: a) SCCs accrete mass from their surroundings and b) that the lifetime of SCCs is proportionally longer than Protostellar clumps at the same mass below 470 M and shorter than Protostellar clumps at masses above 470 M . If mass accretion is the sole explanation then for an average clump free-fall timescale, M ∼ 200 − 440M Myr −1 is required. If the variation in SCC and Protostellar clump lifetimes is the sole explanation, then their ratio of lifetimes should scale as M −0.4 for 100 ≤ M ≤ 1000 M . In reality, both processes likely occur although it is not possible to disentangle them with the current data.