ALMA Observations of Fragmentation, Substructure, and Protostars in High-mass Starless Clump Candidates

The initial physical conditions of high-mass stars and protoclusters remain poorly characterized. To this end, we present the first targeted ALMA Band 6 1.3 mm continuum and spectral line survey toward high-mass starless clump candidates, selecting a sample of 12 of the most massive candidates ( ) within . The joint array maps have a high spatial resolution of ( , θsyn ≈ 0.″8) and have high point-source mass-completeness down to at (or column density sensitivity of ). We discover previously undetected signposts of low-luminosity star formation from and bipolar outflows and other signatures toward 11 out of 12 clumps, showing that current MIR/FIR Galactic plane surveys are incomplete to low- and intermediate-mass protostars ( ), and emphasizing the necessity of high-resolution follow-up. We compare a subset of the observed cores with a suite of radiative transfer models of starless cores. We find a high-mass starless core candidate with a model-derived mass consistent with when integrated over size scales of . Unresolved cores are poorly fit by radiative transfer models of externally heated Plummer density profiles, supporting the interpretation that they are protostellar even without detection of outflows. A high degree of fragmentation with rich substructure is observed toward 10 out of 12 clumps. We extract sources from the maps using a dendrogram to study the characteristic fragmentation length scale. Nearest neighbor separations, when corrected for projection with Monte Carlo random sampling, are consistent with being equal to the clump average thermal Jeans length ( ; i.e., separations equal to ). In the context of previous observations that, on larger scales, see separations consistent with the turbulent Jeans length or the cylindrical thermal Jeans scale ( ), our findings support a hierarchical fragmentation process, where the highest-density regions are not strongly supported against thermal gravitational fragmentation by turbulence or magnetic fields.


INTRODUCTION
High-mass stars (M * > 8 M ) strongly influence the evolution of galaxies and the ISM, yet many fundamental questions remain to be answered concerning the incipient phases of high-mass star formation (e.g., Beuther et al. 2007;Tan et al. 2014;Motte et al. 2017). Observational constraints on the initial physical conditions of protocluster evolution are a nec-the physical properties of representative samples of starless molecular cloud "clumps" 1 .
Recent blind surveys of dust continuum emission at (sub-)millimeter and far-infrared (FIR) wavelengths of the Galactic Plane have identified large statistical samples of clumps, enabling the discovery of those in the earliest evolutionary phases. Such surveys include the Bolocam Galactic Plane Survey 2 (BGPS; Aguirre et al. 2011;Rosolowsky et al. 2010;Ginsburg et al. 2013) at 1.1 mm, ATLASGAL 3 at 870 µm (Schuller et al. 2009;Contreras et al. 2013;Csengeri et al. 2014), JCMT Galactic Plane Survey 4 at 850 µm (Eden et al. 2017, JPS), and Herschel Hi-GAL at 70, 160, 250, 350, and 500 µm (Molinari et al. 2010(Molinari et al. , 2016. Starless clump candidates (SCCs) are identified by cross-matching clump catalogs to catalogs of star formation indicators and selecting clumps unassociated with any indicators. These indicators include 70 µm compact sources, color-selected young stellar objects (YSOs), H 2 O and CH 3 OH masers, and UCHII regions in Svoboda et al. (2016) for the BGPS and in Yuan et al. (2017) for ATLASGAL, in total identifying more than 2 × 10 3 SCCs in the inner-Galaxy. In addition, more than 10 4 clumps without 70 µm sources have been identified from the Hi-GAL survey (Traficante et al. 2015;Elia et al. 2017). In this study we aim to systematically study a representative sample of the highest mass SCCs within 5 kpc in order to understand the fragmentation characteristics at high-spatial resolution, identify potential high-mass starless cores (M 30 M , R 0.1 pc), and search for previously undetected low luminosity protostellar activity.
A variety of physics, including thermal gas pressure, turbulence, magnetic fields, and the geometry of filaments and density gradients, likely play a role in the fragmentation of molecular clouds and the resultant dense core populations. Recent high-resolution observations with millimeter and submillimeter interferometers of high-mass clumps with little sign of star formation reveal significant fragmentation at the early stage of cluster formation (Zhang et al. 2009;Wang et al. 2011Wang et al. , 2012Wang et al. , 2014Zhang et al. 2015;Lu et al. 2015;Beuther et al. 2015a;Sanhueza et al. 2017). These studies found that the most massive fragments in the clumps are at least ten times greater than the thermal Jeans mass, indicating that additional support from turbulence and/or magnetic fields are required. Most of these studies focused on individual clumps and typically have not had the sensitivity to adequately detect fragments of a thermal Jeans mass (detections of 2 M at 4σ rms ). In contrast, the fragmentation scales in nearby molecular clouds have been 1 In this paper we use the term "core" to refer to a dense gas structure that is ∼ 0.1 pc in size and likely to form a single or bound multiple stellar system. Such cores are embedded within larger scale "clumps" that are dense gas structures likely to form a stellar association or cluster, and are ∼ 1 pc in size and ∼ 10 2 − 10 4 M in mass (c.f. Bergin & Tafalla 2007). 2 Data products can be downloaded from https://irsa.ipac.caltech.edu/data/ BOLOCAM GPS/ 3 http://www3.mpifr-bonn.mpg.de/div/atlasgal/ 4 http://apps.canfar.net/storage/list/JPSPR1 studied extensively with recent notable analyses towards Serpens (Friesen et al. 2017), Orion Integral Shaped Filament (Kainulainen et al. 2017), and Perseus (Pokhrel et al. 2018). These studies find support for hierarchical, scale-dependent fragmentation with separations corresponding to a range between thermal Jeans fragmentation and thermal filamentary gravitational fragmentation. It is not understood how these results extend towards earlier evolutionary stages in massive SCCs which are the focus of this work.
Publicly available millimeter and FIR Galactic Plane survey observations do not have sufficient angular resolution at ∼ 20 − 30 (∼ 0.5 pc at 4 kpc) to study the sub-structure and dense core properties in distant SCCs. The high-mass pre-stellar core candidate G028-C1S (M c ∼ 60 M ) studied in Tan et al. (2013) for example was only identified as protostellar until interferometric followup of outflow tracers (Tan et al. 2016;Feng et al. 2016a). High-mass SCCs remain largely unstudied at high-spatial resolution owing to their historical difficulty in identification and typically large heliocentric distances, with only a handful of studies on individual objects to date (Beuther et al. 2015b;Sanhueza et al. 2017). In particular the high-mass starless clump candidate "MM1" of IRDC G28.23-0.19 (Sanhueza et al. 2013) has been studied in detail to determine that it is devoid of star formation indicators, including 3.6 − 70 µm point sources, H 2 O and CH 3 OH masers (Wang et al. 2006;Chambers et al. 2009), and radio continuum (Battersby et al. 2010;Rosero et al. 2016). The global physical properties of G28.23-0.19 MM1 (corresponding to BGPS catalog clump number 4649) are similar to the average properties of the SCCs presented in this work. G28.23-0.19 MM1 is high-mass, cold, compact, and dense (i.e., M cl ≈ 1.5 × 10 3 M , T K ≈ 12 K, R = 0.6 pc, n ≈ 3 × 10 4 cm −3 ; Sanhueza et al. 2017). However the sensitivity and sample size of dense cores are not sufficient for a precise measurement of the fragmentation scale, and represents only a single clump. In this survey we present observations on a sample of 12 clumps that are selected from a blind Galactic Plane survey and are among the most massive SCCs.
Existing large samples of SCCs have been primarily identified through the non-detection of coincident Hi-GAL 70 µm sources (Veneziani et al. 2013;Traficante et al. 2015;Svoboda et al. 2016;Elia et al. 2017) which is less affected than shorter wavelength 8 µm or 24 µm observations by both local extinction and from contamination of evolved stars (principally those on the asympototic giant branch). The completeness of the 70 µm maps to protostar bolometric luminosity, L bol , is affected by the survey depth and complex structure in the foreground and background emission which hinders the clear identification of compact sources. Svoboda et al. (2016) calculate the L bol completeness function for Hi-GAL 70 µm compact sources associated to BGPS clumps and the respective distribution of heliocentric distances and find that for clumps with low 70 µm backgrounds (∼ 500 − 1000 MJy sr −1 ) the 90% completeness limit is L bol = 50 L (see §3.2.4 in Svoboda et al. 2016), which is greater than 95% of YSOs in the Gould's Belt (n.b. median 1 L ; Dunham et al. 2014).
Faint 24 µm sources coincident with the clump column density peaks towards 9/18 of 70 µm dark SCCs suggest likely embedded intermediate-mass star formation that is undetected in the Hi-GAL observations (Traficante et al. 2017). However, it is currently unknown what degree of star formation has initiated, if at all, in SCCs without sensitive and unambiguous tracers of protostellar activity such as bipolar molecular outflows. Systematic observations of SCCs at high-resolution are necessary to determine what degree (if any) of low-luminosity star formation has begun in SCCs, with important implications for the protostellar accretion history.
The principal theories of high-mass star formation in dense Galactic molecular cloud clumps are the monolithic collapse of turbulent cores in virial equilibrium (McKee & Tan 2002, 2003Hosokawa & Omukai 2009) and the accretion of subvirial cores through gravitionally-driven cloud inflow (Smith et al. 2009;Hartmann et al. 2012). The latter replace the competitive Bondi-Hoyle accretion of cores (Bonnell et al. 2001;Wang et al. 2010) with cores fed the gas reservoir through inflowing streams. The turbulent core model predicts monolithic high-mass starless cores whereas the competitive model predicts a fragmentation of cores near the thermal Jeans mass. The existence of high-mass starless cores is a key distinction between these models, yet few, if any, observational candidates are known (Kong et al. 2017), and some promising candidates have revealed embedded protostellar activity upon further observational investigation (i.e. G028- C1S Tan et al. 2013C1S Tan et al. , 2016Feng et al. 2016b). Irrespective of the specific theoretical model, measurements of the fragmentation properties at early evolutionary phases provide valuable observational constraints on the initial physical conditions of high-mass star and cluster formation. To this end, we perform a systematic search for high-mass starless cores towards massive SCCs with the Atacama Large Millimeter/submillimeter Array.
In this paper we present a systematic survey of 12 high-mass SCCs at sub-arcsecond resolution. We present our sample selection, observational setup, and data reduction methodology in Section 2. We describe detections of previously unknown, low-luminosity protostellar activity in Section 3 and the modeling of continuum sources in Section 4. We measure and analyze the fragmentation scale between sources in Section 5, discuss the implications in Section 6, and report our conclusions in Section 7. In Appendix C we include a detailed description of the setup and computation of radiative transfer models analyzed in Section 4.

OBSERVATIONS 2.1. Sample Selection
We have identified SCCs through the combined catalogs and images of primarily two dust continuum Galactic Plane surveys: (1) an evolutionary analysis of BGPS 1.1 mm (Svoboda et al. 2016, hereafter S16), and (2) comparison of the Peretto & Fuller (2009) infrared dark cloud (IRDC) catalog with Hi-GAL images (Traficante et al. 2015). The BGPS observed between −10 • < < 90 • with |b| < 0.5 • (expanding to |b| < 1.5 • for selected ) at λ c = 1.12 mm with a θ hpbw = 33 synthesized angular resolution. In the re-gion 10 • < < 65 • the BGPS has been compared to a diverse set of a observational indicators for star formation activity (S16). These indicators include compact 70 µm sources, mid-IR color-selected YSOs, H 2 O masers, Class II CH 3 OH masers, extended 4.5 µm outflows, and UCHII regions. From the sample of more than 2500 SCCs in the combined samples of S16 and Traficante et al. (2015), we target the 12 highest mass SCCs within d < 5 kpc. Point sources at 70 µm were identified by visual inspection in S16 and by an automated extraction in Traficante et al. (2015). Three clumps (G28565, G29601, and G309120 which were initially determined from the automated extraction to be dark at 70 µm, upon closer scrutiny by visual inspection show association to weak sources. Among the 12 ALMA targets, nine have no detectable point source emission from Hi-GAL 70 µm (flag 0 in S16), two have low-confidence or marginal detections (flag 4, G28565 & G29601), and one has bright, compact detection (flag 1, source G30912). We emphasize that starless clump candidates are designated based on the observational data sets and identification techniques used, and that these factors are reflected in the completeness and purity of the resulting catalogs of SCCs. Table 1 details the target positions and velocities. Table 2 details the physical properties of the sample and Figure 1 shows images of the clumps at wavelengths from 8 µm to 350 µm.
The clump average physical properties in S16 are shown in Figure 2, plotting peak mass surface density Σ cl,pk (at θ hpbw = 33 ) and total mass M cl , for sources with wellconstrained distances less than d < 5 kpc and 10 • < < 65 • . Protostellar clumps and SCCs are plotted, where SCCs have quiescent background emission and no detected compact sources from the Hi-GAL 70 µm images (flag 0, see S16 §3.2.4). Protostellar clumps are typically higher in both mass and mass surface density compared to SCCs. The ALMA targets are shown, occupying the highest Σ cl,pk and M cl portions of the SCC distribution where typical values for the sample are M cl ∼ 800 M and Σ cl,pk ∼ 0.1 g cm −2 (measured over ∼ 0.6 pc scales). G28539 in particular stands out as the most massive clump in the sample at M cl ∼ 3 × 10 3 M and also the highest peak mass surface density. Assuming a star-formation efficiency of sf = 0.3 and a standard stellar initial mass function (IMF; Kroupa 2001), a 320 M clump meets the criteria of forming a 8 M star (see S16 §6.1). All of the clumps that comprise this sample are above this mass threshold, and are similarly above the massradius relationship for high-mass star formation proposed by (Kauffmann & Pillai 2010) 5 . However in practice it is difficult to assess the high-mass star formation potential of clumps beyond these simple heuristics. It should be kept in mind that if the star formation efficiency of the targets is substantially 5 Kauffmann & Pillai (2010) define the prescription M ≤ 580 M (R/pc) 1.33 for a clump to form high-mass stars. The prefactor has been scaled for consistency with the dust opacity used in this work. For radius R = 0.8 pc, this relation yields a mass threshold of M ≈ 430 M . Figure 1. Mid-and far-infrared 3 × 3 maps of the clumps in the survey sample, showing GLIMPSE 8 µm, MIPSGAL 24 µm, and Hi-GAL 70 µm and 350 µm. The ALMA Band 6 single pointings target the peak flux positions derived from the BGPS 1.1 mm observations. The inner and outer red circles show the 50% (27 ) and 20% (40 ) power points of the primary beam for the ALMA 12 m array images. Clumps from Svoboda et al. (2016) were selected to have no detected indicators of star formation activity such as embedded 70 µm sources, H2O and CH3OH masers, and UCHII regions. Clumps from Traficante et al. (2015) were selected to be 70 µm dark using an automated extraction, one of which shows a marginal detection and two of which show clear detections upon visual inspection. Note that G30120 at b ≈ −1.1 • is outside the MIPSGAL survey and does not contain Spitzer 24 µm data. lower than the typical assumed value of sf = 0.3 then they are unlikely to form high-mass stars.

ALMA Band 6
As part of Atacama Large Millimeter/submillimeter Array 6 (ALMA) Cycle 3 program 2015.1.00959.S, we observed 12 clumps in Band 6 in a compact configuration (C36-2; joint 12 + 7 m array baselines range from ∼ 9 − 450 m). Data were taken between 3 − 20 March, 2016, for the 12 m array and between 30 April to 19 August, 2016, for the 7 m array. Including time for calibration and overheads, the 12 m ar-ray observations lasted for approximately 12 hr, with typical precipitable water vapor of 1.5 mm. Titan and J1733-1304 were used as flux calibrators, J1751+0939 to calibrate the bandpass, and J1743-0350 and J1830+0619 to calibrate the time-dependent gains. Identical 1 hr scheduling blocks were configured to interleave and observe all 12 targets within the same block, and because sources are within a 5 • radius on the sky (22. • 7 < < 30. • 9), the same calibrators can be used. Thus due to nearly identical observing conditions, the individual maps have similar uv-coverage, atmospheric noise, and beam size.
Positions for the sample were chosen from the BGPS 1.1 mm continuum peak flux density position, and compared for consistency with the ATLASGAL 870 µm peak emission and Hi-GAL 70 µm peak absorption (when present) posi-R = 0.1 pc 1 pc 10 pc Figure 2. Peak mass surface density Σ cl,pk versus total mass M cl . Values are derived from the BGPS at 1.1 mm (θ hpbw = 33 ) for clumps with well-constrained distances d < 5 kpc. Starless clump candidates (blue points, contours), protostellar clumps (orange contours), and the ALMA sample (blue stars). Total masses of the sample range between M cl ≈ 400 − 3000 M and Σ cl,pk ∼ 0.1 g cm −2 . The dashed lines show Σ cl,pk as a function of M cl for constant radii at 0.1 pc, 1.0 pc, and 10 pc.
tions. The Band 6 receiver was configured in dual-polarization mode with lower and upper sidebands centered near 215 and 230 GHz, respectively. The observations targeted each clump peak with a single pointing with half-power beam width (HPBW) of the measured primary beam 26. 6 (∼ 0.5 pc at d = 4 kpc) and 20%-power beam width of 40 (∼ 0.8 pc), the effective limit of the 12 m array field of view.
2.2.1. ALMA 1.3 mm continuum reduction Data reduction was performed using CASA (version 4.7.134-DEV, r38011, for consistency with QA2 delivered products). Line-free continuum visibilities were created by flagging channels contaminated by spectral lines, where the input spectral windows were further visually inspected to check for emission at unexpected velocity ranges, partitioned out into a new measurement set with the split task, and channel averaged to 25 MHz to avoid bandwidth smearing. Together, this yields ≈ 3.5 GHz of dual polarization continuum bandwidth. The continuum image root mean square (RMS, σ rms = n I 2 n /n) is measured for each CLEANed image within a region that excludes identified emission using the casaviewer tool. None of the images are dynamic range limited with peak image intensity divided by the RMS less than 200. We estimate the fiducial mass sensitivity given T K (NH 3 ) ≈ 12 K, thermally coupled gas and dust (T d = T K ), and OH5 dust opacity κ(λ = 1.3 mm) = 0.899 cm 2 g −1 ). The methods for deriving dust mass values from the continuum emission are discussed in more detail in §4. The joint 12 + 7 m continuum was then iteratively CLEANed with manual masking using the tclean task using the multiscale deconvolver and a robust weighting of 1, down to a brightness threshold of 2 − 3σ rms . An image cell size of 0. 1 was used for all continuum and spectral line maps. Self-calibration was not applied because the brightest sources in the image are only a few mJy, and not sufficiently bright such that a conservative self-cal produces a noticeable improvement without also increasing the image noise. The resultant images have a synthesized beam size of θ maj ≈ 0. 85 by θ min ≈ 0. 75 (0. 8 angular diameter yields 2800 − 3800 au at d = 3.5 − 4.8 kpc). The continuum images are shown in Figure 3.

ALMA spectral line reduction
The flexibility of the ALMA correlator enabled simultaneous observation of several molecular line transitions. Table  3 reports the details of the correlator configuration. We observed nine spectral windows (SPWs) with one wide-band, low-spectral resolution window centered at 233.8 GHz and eight high-spectral resolution windows centered on lines of interest. Table 4 reports the transition quantum numbers, rest frequencies, and upper energy levels (E u /k). The SPWs containing the H 2 CO and CH 3 OH transitions have a spectral resolution of 0.34 km s −1 and the other line SPWs have 0.68 km s −1 resolution. Line rest frequencies were taken from a combination of the SLAIM 7 (F. J. Lovas, private communication, Remijan et al. 2007) and the CDMS (Müller et al. 2005) online spectroscopic databases. The line SPWs from the 12 + 7 m arrays were jointly imaged using the CASA task tclean with a Briggs robust parameter of 1.0, cell size of 0. 1, and re-gridded to common spectral resolutions listed in Table 4. We find a typical RMS noise levels in the image cubes of 1.8 mJy/(km/s) (i.e., 2.2 mJy per 0.68 km s −1 channel or 3.0 mJy per 0.34 km s −1 channel) or 71 mK/(km/s) when converted to brightness temperature units (HPBW beam size of 0. 85 × 0. 75).
In this work we inspect the line image cubes for detection of emission, and for the presence of outflows traced by CO and SiO, but do not CLEAN the data cubes. Due to the lack of full uv-coverage, the CO maps in particular show strong effects of spatial filtering near the systematic velocities that make the deconvolution process complex and error prone. Detailed analysis of the spectral line data is left to a future work. Table  6 lists the detection flags per target for the continuum, molecular lines, and outflows. Features are considered detections if they have peak intensities ≥ 7σ rms ("D"), weak detections if between 5 − 7σ rms ("W"), and non-detections if ≤ 5σ rms . Targets that exhibit bipolar outflows in CO or SiO are flagged "B" (discussed below in §3.2). The dense gas features revealed in the continuum maps clearly show hierarchical structure, with bright ridges, filaments, and cores contained within larger, lower surface brightness features. Given the complexity within the maps and the systematic uncertainties of imaging, we compare the continuum images to an additional measure of gas column density at comparable resolution, MIR extinction. For appropriate configurations of distance and the MIR radiation field, clumps can appear associated with 8 µm absorption features (EMAFs), where high column densities at close distances typically yield the strong MIR shadows that identify infrared dark clouds (IRDCs). MIR extinction mapping has the dual advantages of comparatively high-resolution, insensitivity to dust temperature, and lack of spatial filtering. We use the Spitzer GLIMPSE (Benjamin et al. 2003;Churchwell et al. 2009) IRAC Band 4 (λ c = 7.9 µm, 2 FWHM) mosaic to show the EMAF contrast. Figure 4 presents a map of the flux density S 8 with the ALMA 230 GHz continuum for source G24051 overlaid. The dense gas structures observed in the millimeter continuum show a remarkable consistency when compared with the column density features inferred from the of MIR contrast. This holds similarly true for the other clumps in the sample, as all show at least some MIR extinction. Qualitatively this good correspondence supports the fidelity of the emission structure detected in the ALMA maps. . Comparison between the ALMA 230 GHz continuum (black lines) and IRAC 8 µm intensity S8 (color map, inverted) for clump G24051. Good spatial correspondence is observed between the sub-structure in the ALMA continuum and the highest extinction features in the GLIMPSE 8 µm map (θ fwhm ≈ 2 ). The continuum images are shown without correction for primary beam attenuation for visual display purposes, and the contours are at steps of 2, 3, 5, 10, 20, and 40σrms. The dotted lines show the 50% and 20% power points of the ALMA primary beam.

Core identification and dendrogram
In order to analyze the fragmentation scale we first identify dense gas sub-structures using a segmentation algorithm. The nature of the tree data-structure in the dendrogram algorithm makes it well-suited to identifying and categorizing structure in images with hierarchical structure (see Rosolowsky et al. 2008), as opposed to a simpler segmentation algorithms, such as that done with a seeded watershed algorithm (e.g., CLUMPFIND; Williams et al. 1994). We use the open source Python software library astrodendro to create the dendrogram and catalog of cores. The dendrogram has three principal tunable parameters, defining a minimum threshold value v min setting the floor or outer boundary of each tree, the minimum contrast or step size δ step between nodes, and the minimum area Ω min . Because the noise varies considerably across the primary beam of each image, we apply the dendrogram to maps that have not been corrected for the weight of the primary beam. This effectively works to identify features with outer contours of constant statistical significance across the field of view, rather than outer contours of constant flux. Sources are extracted out to the limit of the maps, set to the 20% power point of the primary beam. We choose conservative values for each parameter, using v min = 3σ rms , δ step = 3σ rms , and Ω min = Ω bm , applied to the unmasked images. Sources (i.e. leaves, or nodes without children) are then sub-selected to meet the criteria that the peak flux is > 5σ rms . In total, we identify 67 sub-structures for the sample of 12 clumps. Figure 5 shows the dendrogram extracted dense gas sub-structures in each clump. Table 5 catalogs the measured positions, sizes, and flux densities of the sub-structures. We find an average number of sub-structures per clump of N src = 5.6 (median 6), with the maximum N src = 11 in G24051 and minimum N src = 1 in G23605. G23605 is the only clump with N src < 3 and is thus not included in the source nearest neighbor distance analysis. Figure 6 presents N src per clump versus Σ cl,pk , where a tentative increasing trend is observed, where clumps that have high Σ cl,pk are more fragmented than lower Σ cl,pk . While in theory the distribution of integrated flux densities can be analyzed to measure a CMF, large observational uncertainties exist in practice that complicate its interpretation. The principal contributor arises from at least a factor of three uncertainty in T d ∼ 6 − 35 K (∼ 10× uncertainty in M ), due to uncertainty in the ISRF, local extinction, and uncertainty in the protostellar activity of each source. From single wavelength observations we do not have enough information to construct spectral energy distributions (SEDs) and measure average line-of-sight dust temperatures. Other significant systematics also arise from uncertainty in the missing flux density due to spatial filtering by the interferometer, dust opacity (δκ/κ ≈ 50%), kinematic-derived heliocentric distance (δd /d ≈ 15%), and the aperture or source boundary used to extract S ν . For these reasons, we shall leave the study of the characteristic fragmentation mass and the CMF in SCCs to a future work utilizing complementary JVLA NH 3 observations that will provide both gas kinetic temperature and kinematic information (Svoboda et al. in prep.). The characteristic fragmentation length scale, on the other hand, can be inferred directly from the distribution of angular separations between sources with assumptions on how to correct for geometric projection.
3. PROTOSTELLAR ACTIVITY In this section we describe new evidence for protostellar activity, and in Sect. 5 we perform an analysis of the fragmentation scale from the sub-structure detected in the continuum. With the improved sensitivity and resolution of ALMA, multiple indicators of protostellar activity are observed for the first time. In particular bipolar molecular outflows detected in CO J = 2 → 1 and SiO J = 5 → 4 provide unambiguous evidence of embedded protostellar activity. The detection of molecular transitions with comparatively high upper excitation temperatures (E-CH 3 OH 4 2,2 → 3 1,2 , E u /k = 45.5 K; p-H 2 CO 3 2,2 → 2 2,1 , E u /k = 68.1 K) and detection of bright, compact continuum emission (unresolved on scales smaller than 3000 au) are also suggestive of embedded, low-L bol protostellar activity. Together, these data provide a clear indication of embedded protostars towards 11 out of 12 clumps.

Compact continuum sources
Numerous high SNR (signal-to-noise ratio, S ν /σ rms ), point-like sources are observed in the continuum images . Dendrogram extracted dense gas sub-structures (orange contour) over-plotted on the ALMA 12 + 7 m array jointly deconvolved 230 GHz line-free continuum images. Elliptical sources are visualized (red ellipses). Sub-structures are labeled by their catalog number from Table 5. The maps are uncorrected for primary beam attenuation for visual display purposes. The color scale ranges from −0.15 to 1.0 mJy beam −1 on a linear scaling. The scalebar (cyan) visualizes 0.3 pc at the clump heliocentric distance. The dashed circle shows the half-power beam width (27 ) and the image extends down to the 20% power point (40 ).
ure 3). We speculate such sources originate from the dense, centrally heated inner-envelopes of embedded protostars. In §4 we investigate whether the compact continuum sources are inconsistent with radiative transfer models of dense, starless cores.
We designate continuum sources as "compact" if they are unresolved or are marginally resolved on the scale of the ALMA synthesized beam θ syn ≈ 0. 8. Continuum sources are determined to be unresolved if a Gaussian fit to the image plane data using the CASA task imfit reports an deconvolved angular sizes θ dec θ syn . The deconvolved Gaussian FWHM are determined from subtracting the synthesized HPBW in quadrature from the fitted width, i.e. θ dec = θ 2 fit − θ 2 syn . These angular widths correspond to physical sizes of 1500 au at heliocentric distances of d ≈ 4 kpc. The brightest compact sources have typical peak flux densities between S 1.3,pk ≈ 1 − 7 mJy beam −1 . All clumps aside from G28539 host a compact source with S 1.3,pk > 1 mJy beam −1 . Indeed, sources G23605 and G30120 host compact sources, even though they show limited fragmentation otherwise. While lacking extended continuum emission, the compact source G23605 S1 in has clear association with emission from multiple molecular species (C 18 O, H 2 CO, CH 3 OH) at the LSR velocity of the clump, determined from the NH 3 emission (32 , 0.7 pc resolution; S16). G30120 S1 is a compact source near the eastern edge of the field with a strong CO outflow and other molecular detections. . Peak clump mass surface density from the BGPS 1.1 mm data versus the number of leaves (i.e. dendrogram leaves) per clump from the ALMA observations. The data hints at an increasing trend of higher mass surface density clumps associated with a higher degree of fragmentation.
pens at 230 GHz with CARMA. The envelope masses range between M env = 0.5 − 20 M (median M env = 3.7 M ) and with integrated flux densities between S 1.3 = 1.4 × 10 1 − 4.0 × 10 3 mJy (median S 1.3 = 120 mJy) and deconvolved size scales between D = 400 − 3000 au (median D ≈ 700 au). With a heliocentric distance of d = 415 ± 25 pc to Serpens (Dzib et al. 2010), the 120 mJy median source flux density and 700 au size measured by Enoch et al. (2011) correspond to 1.2 mJy and 0. 18 when scaled to a fiducial distance of 4 kpc. If there are low-to intermediate-mass Class 0 YSOs with similar physical properties in these SCCs as in Serpens, then they would be consistent with the observed bright ( 20σ rms ) unresolved point continuum sources. This is further supported by the frequent coincidence of outflows towards such sources, discussed in section §3.2. To determine whether the observed compact continuum sources are consistent with starless cores (∼ 0.1 pc) embedded within the mapped clumps (∼ 1 pc), in §4 we compare a subset of the observations to radiative transfer models of starless cores. Some continuum sources without molecular line detections may be background galaxies. Deep surveys performed with ALMA (Hatsukade et al. 2013;Carniani et al. 2015) have determined source counts of background galaxies at 1.3 mm. The number of sources expected in the images with flux densities greater than N (S 1.3 > 0.3 mJy) 3 over the 12 fields, measured as HPBW area for each pointing, outside of which the degraded sensitivity yields negligible background sources. This represents approximately 5% of the detected sources, and thus does not have a significant effect on the calcula-tion of the nearest neighbor separations or other estimated distributions of core properties.

CO & SiO Outflows
Ordered, bipolar molecular outflows driven by protostellar accretion provide a sensitive and unambiguous detection of embedded protostellar activity (see reviews by Arce et al. 2007;Frank et al. 2014). In this section we describe the properties of outflows detected with ALMA in CO J = 2 → 1 and SiO J = 5 → 4.
Outflows are identified through visual inspection of the CO datacubes in conjunction with the 1.3 mm maps overplotted. While the emission structures in the CO cubes are complex, bipolar outflows are clearly apparent as paired linear emission structures. These features are identified as linear features radiating from the same location with highly ordered red and blue velocity components that are detected over many velocity channels ( 10 km s −1 , 15 channels). Outflow candidate features with only a single red or blue component are also observed, but due to the greater ambiguity in identification these are not regarded as clear signatures of star formation activity. The CO outflows are generally highly ordered in position and velocity, but spatial filtering of bright, extended emission and self-absorption near the source systemic velocity complicate the identification of low-velocity (|v| 1.5 km s −1 ) outflow components. Higher velocity components of the spectra also suffer both self-absorption from foreground CO clouds and confusion with bright Galactic emission, which can bias measurements of the maximum outflow velocity to lower values. Analysis of an example outflow in G24051 is presented in Appendix A.
We find that 9 out of 12 clumps are associated with bipolar CO outflows and 16 outflows in total are observed. We also find that 3 out of 12 clumps are associated with bipolar SiO outflows and 4 outflows in total are observed. The clumps with outflows are reported in Table 6. Pairs of CO outflows originating from the same continuum source are also observed, as seen in G23297 S2 and G29601 S1, which point to unresolved protostellar multiple systems. Figure 7 presents the ALMA joint 12 + 7 m array CO J = 2 → 1 integrated intensity maps for blue-and red-shifted velocity components.
We also detect SiO emission towards several more continuum sources and positions without clear signs of ordered bipolar outflows. SiO emission detection is a strong indicator of protostellar activity because of its origin in high-velocity shocks driven by protostellar outflows (Schilke et al. 1997). However, recent work has shown that low-velocity shocks ( 10 km s −1 ) created by colliding flows may produce substantial distributed SiO emission (Jiménez-Serra et al. 2010;Nguyen-Lu'o'ng et al. 2013;Louvet et al. 2016). Thus, considered by itself, a detection of relatively narrow linewidth (∆v 10 km s −1 ) SiO J = 5 → 4 emission is not an unambiguous indicator of star formation activity. Maps of SiO integrated intensities are presented in Appendix B.
The 70 µm dark clump G28539 (upper-left corner of Fig. 3 & 7) shows no clear sign of CO or SiO outflows, and thus remains a starless clump candidate at the improved sensitivity of ALMA. Several indirect tracers of star formation are observed towards G28539 however and we discuss these in turn.
Moderately high-excitation molecular lines (E u /K 50 K) are unlikely to be excited in the cold 10 K gas expected to be found in starless cores and quiescent clump gas. Detection of such lines in our observations are thus indirect evidence of embedded protostars, although as discussed in §3.2 it is possible some of these lines are excited from low-velocity shocks originating from colliding flows. In G28539, a compact source of weak emission CH 3 OH and H 2 CO 3 2,2 → 2 2,1 is detected. These features are not coincident with continuum emission and may originate from non-protostellar shocks, shocks of undetected protostellar outflows, or of embedded protostellar cores that are below our detection limit. Similarly, a compact source of SiO is also detected does not coincide with any continuum emission feature (it may be seen on the west side of the field in Fig. 15).
There exist weak 24 µm sources in the vicinity of the clump boundaries as defined by the 350 µm and 500 µm emission, notably a faint source within the extinction feature ∼ 1 east of the ALMA field (see 24 µm panel in Fig. 1), a brighter source on the south-eastern outskirt of the clump, and a marginal feature coincident with the continuum source in the NW edge of the ALMA field of view. Because of the substantial con-tamination from evolved stars, 24 µm emission alone is not a robust indicator of protostellar activity. If these sources are indeed protostars associated with the clump then they would be evidence that star formation has begun in G28539.
Deep radio continuum observations when available also provide a diagnostic of star formation activity because they are sensitive to the ionized gas in ultra-and hyper-compact HII regions, ionized winds, and jets from low-to intermediatemass protostars. Rosero et al. (2016) carried out deep Jansky Very Large Array (JVLA) C & K-band observations towards a sample of high-mass clumps which contains source G28539 in the field "G28.53-00.25". The HPBW of the primary beam for the JVLA at C-band is 9.2 at 4.2 GHz (LSB) and 4.2 at 7.4 GHz (USB), with synthesized HPBW resolution of approximately ∼ 0.4 in the A-configuration. Using the radio-continuum to bolometric luminosity scaling relations for protostars in Shirley et al. (2007) (Eq. 3), the measured σ rms = 3 µJy beam −1 sensitivity at d = 4.7 kpc can be converted to a 1σ bolometric luminosity sensitivity of ∼ 30 L , that is reasonably comparable to the PACS 70 µm sensitivity from Hi-GAL. Here a faint point source is detected near the center of the ALMA pointing, detected in both sidebands at moderate significance (8 and 5σ in LSB, USB respectively). The measured in-band spectral index (S ∝ ν +α ) α = −0.65 ± 0.46 favors a non-thermal synchrotron dominated source, but the weak constraint is consistent with thermal free-free emission α = −0.1 at 1.2σ rms . The location 18 • 44 22. 621 −4 h 02 m 00. s 380 (J2000) is not coincident with millimeter continuum or spectral line emission in the ALMA data. Given the lack of a clear association, we conclude that this radio continuum source is likely an extra-galactic contaminant and not an indicator of protostellar activity.
In summary, indirect evidence for star formation exists from two different tracers: (1) 24 µm sources at the edge or outside of the ALMA field of view, and (2) ALMA detections of CH 3 OH and SiO that are not clearly associated with continuum sources. G28539 is the most massive clump in the sample (M cl ≈ 3600 M ) and shows fairly limited signs of fragmentation. After the ALMA observations G28539 is the only starless clump candidate remaining in our sample. It is thus a target of great interest for studying the initial conditions of high-mass star formation.

MODELING CONTINUUM SOURCES 4.1. Starless core models
A diverse range of continuum sub-structures are found to be present in SCCs, from unresolved compact sources, filaments, to lower surface brightness extended emission. In this section we analyze whether cores with bright, unresolved continuum emission on scales < 1500 au (∼ θ syn /2) are necessarily protostellar even without detections of outflows or strong highexcitation molecular lines. We also model whether low-to intermediate-mass starless cores are accurately recovered in the observations and perform detailed modeling of high-mass starless core candidates in clump G28539.
To characterize the continuum features in our images we apply the radiative transfer code RADMC-3D (Dullemond et al. 2012) to self-consistently calculate the equilibrium dust temperature distributions of externally heated starless cores and to produce synthetic images. We follow a similar approach to modeling starless cores as found in Shirley et al. (2005) and Lippok et al. (2016). We apply conventional assumptions for the dust properties (Ossenkopf & Henning 1994;Weingartner & Draine 2001;Young & Evans 2005) and interstellar radiation field (ISRF, Draine 1978;Black 1994). A detailed description of the computed models may be found in Appendix C.
We apply a spherically symmetric Plummer-like function to parametrize the model radial density profile (Plummer 1911;Whitworth & Ward-Thompson 2001;Lippok et al. 2016). The gas density profile n H can be expressed as: (1) for radius r, inner gas density n in , outer gas density n out , flat radius R flat , and power law exponent η (n.b. an isothermal Bonnor-Ebert sphere may be approximate with η = 2; Ebert 1955;Bonnor 1956). The strength of the interstellar radiation field (ISRF) is varied from the local value by a multiplicative scale factor s isrf . We compute 10 4 models, randomly sampling the parameter space by drawing values from a uniform distribution in log-space within the ranges for the parameters n in = 1 × 10 4 − 1 × 10 7 cm −3 , n out = 1 × 10 1 − 1 × 10 3 cm −3 , R flat = 1 × 10 3 − 2 × 10 4 au, and s isrf = 1 − 100, while η = 2.5 − 5.5 is drawn uniformly in linear space. Models are evaluated on a logarithmic radial grid from 2.5 × 10 2 au to 6.0 × 10 4 au. These values are chosen to cover the range of values from the sample of low-and intermediate mass cores in Lippok et al. (2016) but extended to higher n in and smaller R flat . After computing the radiative transfer, the models are ray-traced by RADMC-3D and projected to a fiducial distance of d = 4 kpc.

Model recovery
We find that 53% of the computed models (5268/10 4 ) meet the detection threshold of S 1.3mm,pk > 5σ rms when convolved with a θ = 0. 8 Gaussian beam. The cut in peak flux density has no effect on the recovered distributions of η and n out , and minimal effects on R flat and s isrf , with an increase in the median values by a factor of 1.5 and 1.2, respectively, over the distribution of model cores.
It is important to keep in mind that the suite of model cores is constructed to span the parameter space of relevant values, not to represent an observed or predicted core mass function. We do not use the fractions of detectable cores to infer completeness, but to show the expected range of physical parameters for which cores can be recovered. To estimate this, we sort the models by M and s isrf , counting both the fraction of detectable cores, and regions of parameter space with at least one detectable model (Fig. 8). Computing the detection fraction in this way has the effect of marginalizing over our uncertainty in R flat , η, n in , and n out , that are poorly constrained with our single wavelength maps. Figure 8 (left) shows that at 50% completeness, the cores at s isrf ∼ 3 are recovered for M 4 M , and that this extends down to M ∼ 1 M for the extreme value s isrf = 100. Figure 8 (right) shows that it is possible however to recover lower-mass cores if the ranges of models is restricted to those that are the most compact (where R flat < 3 × 10 3 au, η > 4.5) and have high central densities (n in > 1 × 10 5 cm −3 ). For these compact sources, M ∼ 1 M models may be recovered at s isrf ∼ 3 and down to M ∼ 0.2 M for s isrf = 100.
From these models we can infer that the completeness expected from our point-source sensitivity, ∼ 0.3 M at 6σ rms , is an underestimate if the majority of cores are resolved (see also Appendix A in Beuther et al. 2018). Low-mass cores with extended profiles will thus go undetected with a criteria based on peak intensity, leading to our seemingly shallow limit of ∼ 1 − 4 M . Observationally, we must approach extended emission at low-SNR with caution because there is extended structure in the maps on scales larger than the 12 m primary beam that cannot be adequately CLEANed. For this reason we do not attempt to identify or catalog sources down to the limits of statistical significance for extended and spatially-integrated flux densities but maintain a conservative detection limit based on source peak flux density. The typical integrated flux density of a source is S 1.1 ∼ 1 − 10 mJy (M ∼ 1 − 10 M assuming T d = 12 K and d = 4 kpc), and generally consistent with the thermal Jeans mass M j,th for a uniform medium at the density of the clump, M j,th ∼ 2 M where M j,th ≡ (4π/3)(λ j,th /2) 3 ρ 0 for thermal Jeans length λ j,th and average density ρ 0 (McKee & Ostriker 2007; see §5 for an analysis of the Jeans length λ j,th ).

Synthetic observations with CASA
We now investigate whether the models of starless cores provide adequate fits to the brightness profiles present in the SCCs of this survey. We find that compact sources of continuum emission that are unresolved (i.e. deconvolved sizes 1500 au, ≈ θ syn /2) are poorly fit by models of starless cores. Without multiple wavelength observations or gas kinetic temperature information, the radial dust temperature profiles of the cores are poorly constrained. Because of the substantial systematic uncertainties presented in single wavelength observations and potentially undetected embedded protostars, we do not perform a fit to every continuum source, but select a few characteristic examples for quantitative comparison. We create synthetic observations from the models using the CASA sm module by predicting onto the observed visibilities (gridded beforehand for computational efficiency) and imaged without noise using the same tclean configuration as the observations. This does not introduce a significant effect on the models however because nearly all the flux is concentrated on radii r < 2 × 10 4 au or angular diameters of 8 , appreciably less than half of the 12 m array 27 HPBW, and substantially less than the maximum recoverable scale of 33 from the 7 m array. A subset of models were further tested for consistency, because the aim of this comparison is for an understanding of a few representative sources and not detailed parameter estimation, we do not image the full suite of models; rather, we convolve the models with the angular size of the synthesized beam (θ syn = 0. 8) and convert to radial brightness profiles.

Comparison to observations
We compare the observations and models using a method based on the χ 2 -statistic, where the reduced χ 2 r may be expressed as for degrees of freedom ν, independent measurements i, measurements o i , model values m i , and variances σ 2 . We discriminate between models based on the goodness of fit metric ∆χ 2 r ≡ χ 2 r − χ 2 r,best from Robitaille et al. (2007) and Robitaille (2017). Robitaille et al. apply the heuristic that models with ∆χ 2 r < 3 are considered good fits and rejected as poor fits otherwise. Robitaille et al. further note that the Bayesian likelihood under the assumption of normal errors (i.e., P (D|θ j , M ) ∝ exp −χ 2 /2 for data D, parameters θ j , and model M ) yields too stringent a definition of probability given systematic sources of error in the measurements and poor physical correspondence of the model to nature. This ultimately provides a more conservative criteria for rejecting poor fits as the ∆χ 2 r heuristic likely overestimates uncertainties.
We consider two example starless core candidates, G28539 S2 and S4 (see Fig. 5 and Table 5, because they (1) lack unresolved continuum emission at their center, (2) host no outflows or other indicators of star formation activity, and (3) are relatively isolated such that radial brightness profiles can be adequately extracted. G28539 S2 and S4 are also of interest because they are among the brightest such sources, and thus are good high-mass starless core candidates.
We extract radial brightness profiles for the cores by extracting the integrated flux density within 0. 2 diameter annuli about the central position. Uncertainties in the integrated flux densities are calculated as the δS ν = σ rms Ω ann /Ω bm for the solid angle of the annulus Ω ann and the synthesized beam solid angle Ω bm . The radial brightness profiles of the models are then compared by Equation 4.4 for degrees of freedom ν = r max /0. 2−5 ≈ 15 (maximum radius r max = 3. 5−4. 0). Well-fit models are then selected where ∆χ 2 r < 3. Figure 9 shows the best fit models compared to the observations, and Figure 10 shows the radial brightness profiles with the range of fits. We find that the extended brightness profiles are well fit by the starless core models (χ 2 r,best = 1.1 and 0.12 for S2 and S4 respectively). If the range of models are limited to those that resulted in T d (r = 1 × 10 3 au) = 7 − 13 K in order to be broadly consistent with the clump-average temperature derived from the Hi-GAL SED and GBT NH 3 fits (see also the detailed considerations in Tan et al. 2013), then M S2 = 29 52 15 M and M S4 = 14 34 6.0 M , for the median, maximum, and minimum model mass. With a core star-formation efficiency of 30% it is possible that these cores may form high-mass stars (M * > 8 M ). Assuming a 50% formation efficiency from models regulated by outflows (Zhang et al.  Figure 8. Left: Fraction of the computed models with peak flux densities S 1.3mm,pk meeting the source detection criteria > 5σrms in the images as a function of M (r < 2 × 10 4 au) and s isrf . Lower mass cores meet the criteria for larger values of s isrf . The 50% detection threshold for all models (black line) and 50% detection threshold for models with nin > 10 5 cm −3 (gray line) are shown. Cores in this range with M ∼ 1 − 6 lie above this threshold, depending on s isrf , and are relatively insensitive to the choice in model parameters. Note that for log 10 (s isrf ) ∼ 0.5 (or s isrf ∼ 3) the distribution above M 4 M meet the detection criteria. Right: Detection criteria for "compact" models (R flat < 3 × 10 3 au, η > 4.5) with high central densities (nin > 10 5 cm −3 ), for cases where any model meets the criteria (red) and cases where none do (blue). Compact starless core models with M ∼ 1 M are detectable at s isrf ∼ 3. Regions that are not sampled by the compact subset of models are shown in gray, note that because the maximum nin = 10 7 cm −3 , models with M 10 M are more extended.
2014) the maximum expected stellar mass for S2 could be M * ≈ 26 M . Given the fact that these cores are not associated with outflows in the ALMA data or other high-excitation molecular lines, they are excellent candidates for high-mass starless cores. G29558 S1 represents the class of compact continuum sources in our data set. Analysis of this source is then a test of whether the compact sources are well described by starless core models, or alternatively, likely to host embedded protostars. This continuum source has some surrounding extended continuum emission, does not show clear outflows traced by CO or SiO, but is associated with weak CH 3 OH and p-H 2 CO emission. It is bright with peak flux density 6.6 mJy and is similar to other continuum sources with associated outflows. We find that the models poorly fit the observations, with χ 2 r,best = 23.9 and no models for ∆χ 2 r < 9. The properties are pushed to the extremes of parameter space: s isrf ∼ 100, η ∼ 5.5, and n in 1 × 10 7 cm −3 . The moderate R flat ∼ 5 × 10 3 au is a compromise between the compact and extended components of the brightness profile. The poor model fits to G29558 S1 do not strictly require that it or any other individual source is protostellar (models with n in 10 8 cm −3 and R flat < 10 3 au would likely fit the observations). However, such extreme starless cores are unlikely to be observed in significant numbers in our sample, where ∼ 40% of fragments are compact continuum sources. The free-fall timescale of a core with n in = 10 8 cm −3 would be t ff ≈ 3 × 10 3 yr, and for n in = 10 7 cm −3 would be t ff ≈ 1 × 10 4 yr. These are shorter than the inferred ages from the extent and velocity of the observed outflows, although these have substantial uncertainties. Together, the observed properties of these compact continuum sources are more favorably explained as embedded low-to intermediatemass YSOs, which at ∼ 4 kpc would be both of comparable brightness and unresolved (see §3.1). A detailed analysis of the starless and protostellar core properties and dynamics will follow in a future work incorporating NH 3 data from the VLA observations.

FRAGMENTATION SCALE 5.1. Nearest neighbor separations and Monte Carlo simulations
We characterize the linear fragmentation scale in terms of the nearest neighbor separation δ nns between dendrogram leaves in each clump. Geometric projection of sources in the plane of the sky will systematically decrease δ nns from the true value, δ nns . In this work we employ Monte Carlo random sampling to de-project δ nns statistically. Thus while the uncertainty in δ nns may make constraints for any individual pair of sources quite weak, with prior assumptions on the relative positions of sources, the posterior distribution from the ensemble of all δ nns measurements in our sample of SCCs can be readily constrained.
Monte Carlo sampling is used to draw realizations of relative line-of-sight distances z, computing δ nns for each source from the Cartesian coordinates (x, y, z). We use the hierarchical classification of sources in the dendrogram to discriminate between two methods of drawing z values: (i) iso-
lated sources and (ii) sources with common surrounding emission. If sources are isolated (Case i), forming a tree with a single branch, then for each trial we draw line-of-sight distances from a Gaussian distribution with standard deviation σ z = 0.15 pc (FWHM 0.35 pc), chosen such that the doublesided 2σ z interval is 0.6 pc, which is the approximate diameter inferred from the 8 µm maps (cf. Fig. 1 & 4). If sources are associated within the same branch of the dendrogram (Case ii; ie. they are within a common base iso-contour of emission) then we assume that those sources are connected in a filamentary gas structure with unknown inclination with respect to the observer. For each trial, we draw a common inclination φ for the group, pivoting along the major axis, with the pivot axis fixed to z = 0 at the projected geometric center.
Inclinations are drawn such that the length between the two components with the maximum separation δ max is less than D = 0.6 pc, thus where φ is drawn uniformly within the interval (− arccos(δ max /D), + arccos(δ max /D)). If δ max > D, then φ is drawn uniformly within (−65 • , 65 • ), such that δ nns 1.4 pc to extend out to a typical clump effective radius of R ≈ 0.7 pc. In total there are 17 (26%) isolated sources and 49 (74%) grouped sources. Without more detailed knowledge available, informed from either additional observational data or theoretical simulations, we consider this scheme a conservative way to correct the data for geometric projection. While the assumptions in the correction are simple and imperfect, for brevity we refer to the distributions of MC trials as "projection-corrected" below to distinguish it from the projected data. Extensions of this method may opt to use more sophisticated schemes to group sources beyond common millimeter continuum emission, such as grouping sources through a lower density kinematic tracer or a source-density based clustering algorithm.
With no correction applied, the distribution of projected separations has a median value of µ 1/2 (δ nns ) = 0.083 pc with (16, 84) percentile interval of (0.051, 0.140) pc. To calculate the projection-corrected separations, we compute 1 × 10 4 realizations for each clump, and find µ 1/2 (δ nns ) = 0.118 pc with µ 1/2 (δ nns )/µ 1/2 (δ nns ) = 1.42 and a (16, 84) percentile interval of (0.065, 0.232) pc. For comparison, if we assume that all sources are uniformly distributed within a spherical volume of radius R s the following projection correction may be applied: as is done in Myers (2017). If we assume R s = 0.38 pc from the radius of the 20%-power point of the ALMA primary beam at 4 kpc, then this correction factor would be δ nns /δ nns ≈ 1.84 and δ nns = 0.153 pc, which is larger than the median value computed above from the MC trials by 29%.

Jeans length comparison
To consistently compare δ nns values between clumps with different physical conditions, we scale the values by the clump average thermal Jeans length, the minimum wavelength for gravitational fragmentation in an isothermal, uniform medium.

G28539 S2
G28539 S4 G29558 S1 Figure 10. Example ALMA observed sources fit with the suite of starless core models. The observed radial brightness profiles (black) and the image 1σrms (gray region) are shown with the best fit model (red dashed) and envelope of all models that satisfy χ 2 r − χ 2 r,best < 3 (red dotted). The error envelope is calculated as the ±1σ uncertainty of the integrated intensity within the annular aperture at the angular radius θ. The profiles are truncated to where the source is mostly symmetric. The map RMS is visualized (gray dashed line). Top: χ 2 r,best = 1.1. Middle: χ 2 r,best = 0.12. Bottom: χ 2 r,best = 23.9 (magenta dashed), no models for χ 2 r − χ 2 r,best < 9. The resolved sources G28539 S2 and S4 are well fit by starless core models, while the models fail to fit the high-SNR, unresolved inner component in G29558 S1.
The thermal Jeans length λ j,th can be expressed as (McKee & Ostriker 2007): where c s = kT /µm p is the isothermal sound speed (0.21 km s −1 for T d = 12 K), G is the gravitational constant, and ρ 0 is the average volume density. For the accurate propagation of uncertainties in the calculation of λ j,th , we perform MC random sampling of the relevant observational uncertainties in ρ 0 from the dust mass surface density (ρ 0 = 3Σ/4R) and heliocentric distance. The total (i.e., gas) mass surface density are calculated with for source integrated flux density S ν,int , source solid angle Ω, Planck function B ν (T d ) evaluated at dust temperature T d , opacity per mass of dust κ (λ = 1.3 mm) = 0.90 cm 2 g −1 Ossenkopf & Henning (1994), mean molecular weight µ = 2.33, and dust-to-gas mass ratio f d ≡ (m d /m g ) = 1/110 (values are further described in Appendix C). The fragmentation measured within the ALMA maps is most sensitive within the HPBW (27 ) of the primary beam, so an estimate of ρ 0 within this volume we consider to be the most representative density for the computation of λ j,th . Clump average densities on angular scales (∼ 1 − 2 ) larger than the HPBW likely underestimate ρ 0 . Likewise, imageintegrated flux densities from the 12 + 7 m array data possibly underestimate the Σ from spatial filtering. Due to the unfavorable match in resolution compared to the Hi-GAL 500 µm (θ hpbw ≈ 35 ) or BGPS 1.1 mm (θ hpbw ≈ 33 ), we extract flux densities from the ATLASGAL 870 µm maps (θ hpbw ≈ 19 ) at the position of the ALMA pointing for each clump within a beam-sized 27 diameter circular aperture to measure Σ cl . Use of the single millimeter flux mitigates one systematic uncertainty in choosing between Hi-GAL SED fits with or without the 160 µm band included, or using Hi-GAL SED fits that are over the emission for the full clump rather than the peak at consistent angular resolution. The clump average dust temperatures from SED fits to the Hi-GAL data range from T d = 10 − 14 K, but some systematic uncertainty exists with averaging over larger volumes than the ALMA field of view and choices in including the 160 µm band. We choose a conservative dust temperature distribution by assuming a Gaussian dust temperature distribution T d = 12 ± 2 K (1σ interval). For consistency this temperature is also used for the gas kinetic temperature in c s . We propagate the uncertainty in heliocentric distance based on the distance probability density function (DPDF) from Ellsworth-Bowers et al. (2015) for each clump. All sources are well-resolved to the near kinematic distance, and have a δd /d ≈ 0.15 fractional uncertainty. We sample the distributions for S 870 , T d , and d for each MC trial of ρ 0 in the calculation of λ j,th to combine with a trial of δ nns to compute the quotient δ nns / λ j,th for each clump. The computed median volume densities for the clumps in the sample range between n(H 2 ) = (2 − 6) × 10 4 cm −3 and with associated values of the thermal Jeans length between λ j,th = 0.10 − 0.17 pc (2.1 − 3.5 × 10 4 au). The median of samples from all clumps is λ j,th = 0.135 pc (2.77 × 10 4 au).
No correlation is observed between δ nns / λ j,th and the number of cores/leaves in each clump (Fig. 11).
Probability density functions (PDFs) of δ nns / λ j,th are computed for each clump by performing Monte Carlo ran- dom sampling of the observational uncertainties in λ j,th as described above and sampling the de-projected source separations (see §5.1). Figure 12 (left) shows the distributions of δ nns / λ j,th for each clump sorted in descending order by the number of continuum sources. The separation distributions show a bi-modal tendency with peaks at δ nns / λ j,th ∼ 0.3 and δ nns / λ j,th ∼ 1, and with longtails extending to high values 1.5. The distinct peaks at small values of δ nns / λ j,th (all well-resolved) likely result from closely spaced, connected sources where δ nns is not strongly effected from sampling the inclination distribution. Median values of the distributions range between δ nns / λ j,th = 0.4 − 1.5. The values are generally consistent with the thermal Jeans length, but the high frequency of sources with sub-Jeans separations may indicate hierarchical fragmentation at multiple scales. With the initial fragmentation on the clump scale, a further fragmentation on the "core scale" would proceed on sizes 2 × 10 4 au and densities 3 × 10 5 cm −3 . If such hierarchical fragmentation proceeds principally with two resultant fragments on the core scale, then the second nearest neighbor distance would measure the above level in the hierarchy and recover the spacing of the clump scale. This is supported by a plot of the second nearest neighbor distance δ (2) nns distributions, shown in Figure 12 (right), that shows clumps with more uni-modal distributions, with modes and median values at or slightly above the thermal Jeans length. Median values of the δ (2) nns distributions are greater than those for δ nns but generally fall within a similar range between δ (2) nns / λ j,th = 0.75 − 1.7. We compute PDFs for each clump (see above) and the ensemble distribution composed of all separation measurements from each clump aggregated together (Fig. 13). The ensemble separation distribution is used to define a representative fragmentation scale from the SCCs in this survey. As these clumps are at similar distances and blindly selected from Galactic Plane dust continuum surveys, the measured ensemble sample properties may be used to cautiously infer the properties of the Galactic high-mass SCC population (M cl 10 3 M ). Additional observations are required to directly constrain the properties of SCCs with M cl 10 4 M (if they exist outside of the Central Molecular Zone) or SCCs below the mass range of this sample, M cl 400 M . Figure 13 shows the cumulative distribution function (CDF) for the ensemble of δ nns / λ j,th measurements as drawn from the MC sampling for the projected separations, projection-corrected separations, and relevant scales such as the resolution and primary beam HPBW. The projection-corrected ensemble distribution has a median value of δ nns / λ j,th = 0.82 with a (25, 75) percentile interval of 0.52 − 1.25. The percentiles for δ nns / λ j,th = 0.5, 1, 2, and 3 are, respectively, 23. 6, 63.3, 90.3, and 97.4. Overall, the sample of SCCs show a fragmentation scale that is well characterized by the thermal Jeans length.
A relatively small fraction of the separation distribution is inconsistent with the thermal Jeans length, < 10% for > 2 × λ j,th . The large separations do not result from a single or small number of clumps with consistently large separations, but from isolated individual sources within clumps that show fragmentation near the thermal Jeans length. G30660 and G30912, for example, have a significant proportion of the distribution at large separations (see Fig. 12 left), but do not have peculiar dust temperatures, between T d = 11 − 12 K from Hi-GAL SED fits. This portion of the separation distribution may indicate an additional scale for hierarchical fragmentation where a source of non-thermal support prevents fragmentation at the thermal Jeans scale.
The Jeans length can further take into account sources of non-thermal support, such as turbulence or magnetic fields, by using an effective sound speed c s,eff = c 2 s + σ 2 nt 1/2 (6) through the contribution of a non-thermal velocity dispersion σ nt . From S16, 9 out of 12 clumps have T K measured from NH 3 (at 32 resolution). The measured velocity dispersions (i.e. c s,eff ) determined from the spectral line model fit range between σ(NH 3 ) = 0.50 − 0.95 km s −1 with a median value of 0.65 km s −1 , corresponding to σ nt ≈ 0.62 km s −1 for c s = 0.21 km s −1 at T K = 12 K (where T K = 11 − 14 K).
Replacing c s with c s,eff in Equation 5.2 yields the effective Jeans length, or when turbulence is the dominant source of non-thermal support, the turbulent Jeans length λ j,tu . Because λ j ∝ c s , the increase c s,eff /c s ∼ 2.4 − 4.5 (median 3.1) yields a similar scaling for λ j,tu /λ j,th . In comparison, δ nns / λ j,th = 3 (δ nns / λ j,tu ≈ 0.32) occurs at the 97.4 percentile, and thus while such separations are not absent from the data, they are also not representative of the fragmentation measured within the ALMA maps. The length scale distribution is incomplete beyond δ nns / λ j,th > 3.1, where 10% of the MC trials would have 3D separations greater than or equal to the FOV (40 ).

DISCUSSION
The physical processes regulating fragmentation in molecular clouds remain an open problem in star formation. How much are SCCs supported against gravitationally induced fragmentation from non-thermal forms of pressure, such as magnetic fields (B-fields) and/or turbulence? Individual SCCs have been studied at high-resolution (Beuther et al. 2015b;Sanhueza et al. 2017), but we shall discuss a systematic set of observations on a representative sample of high-mass SCCs. Here we describe the fragmentation characteristics of SCCs in the context of theoretical models of star and cluster formation and compare to existing high-resolution observations of clumps and IRDCs.

Cylindrical Fragmentation in SCCs
As shown in §5, we find that clumps fragment at scales consistent with the thermal Jeans length in SCCs. It is known however that geometry and non-thermal support affect the predicted fragmentation scale, producing deviations from that expected for an isothermal, uniform medium. In this section we discuss how the fragmentation scale observed with ALMA compares to different characteristic length scales.
Filaments are ubiquitous in both observed molecular clouds and simulations (e.g., Barnard 1907;André et al. 2014;Smith et al. 2016), and thus cylindrical geometry is of special significance to dense molecular regions. On larger spatial scales observable in MIR extinction, it is clear that the clump peaks are embedded in filamentary gas structures (see Fig. 1 and G23297 for a good example). An infinite, self-gravitating cylinder is unstable to axisymmetric perturbations or "sausage" instability, where the cylinder fragments at the scale of the fastest growing mode of the fluid instability (Chandrasekhar & Fermi 1953;Ostriker 1964;Larson 1985;Nagasawa 1987;Inutsuka & Miyama 1992). For a pressure confined isothermal gas cylinder of radius R and scale height H = c s (4πGρ c ) −1/2 (where ρ c is the central density of the cylinder) then the fastest growing mode depends on the ratio of R and H (Nagasawa 1987). In the case where R H then λ cyl ≈ 10.8R; and, alternatively where R H, then λ cyl ≈ 22.4H (Nagasawa 1987;Jackson et al. 2010).
Is the isothermal, cylindrical fragmentation scale representative in SCCs? The approximation of SCCs as isothermal is imperfect due to shielding that decreases the temperatures of inner regions, but the assumption is generally more valid than clumps with active HMSF and substantial internal protostellar heating and feedback. Observed aspect ratios of 5 over the full clump extent support the approximation of an infinite cylindrical geometry. The typical radial extent of the SCCs as observed in the MIR extinction maps suggests R ∼ 0.4 pc. Assuming that the cylinder central density is equal to the observed clump peak density (i.e. ρ c = ρ 0 ≈ 3 × 10 4 cm −3 ) then H ∼ 0.02 pc, and thus R/H ∼ 20 roughly satisfies the condition R H. Note that for ρ c equal to the clump peak density then this simplifies to λ cyl /λ j,th ≈ 3.50, or for the median λ j,th = 0.137 pc, λ cyl = 0.480 pc. We find that λ cyl is not representative of the δ nns distribution in SCCs, with the observed δ nns / λ j,th ∼ 1.
Observational studies carried out on larger spatial scales than this work support λ cyl as a characteristic scale in filaments (Beuther et al. 2015a;Friesen et al. 2016). While these studies did not have sufficient spatial resolution to adequately resolve the thermal Jeans length, they probe separations on the clump scale and larger than ∼ 1 pc, as observed with ALMA. This work complements the larger-scale studies by identifying fragmentation on the clump Jeans length at an early evolutionary phase. This is supported by the results of Kainulainen et al. (2013) who with Spitzer MIR extinction mapping find that the molecular filament G11.11-0.12 is well described by filament fragmentation and turbulent λ cyl on δ 0.5 pc and λ j,th on smaller scales. Beuther et al. (2015a) in an analysis of the fragmentation in the star-forming filament IRDC 18223 find a mean fragment separation of δ = 0.40 ± 0.18 pc, consistent with a thermal λ cyl = 0.44 pc of the filament, approximately twice that of λ j,th = 0.07 − 0.23 pc, however the authors note that measures of δ should be considered an upper limit due to the sensitivity and resolution of the data. Friesen et al. (2016) in a survey of the entire Serpens South molecular cloud (as part of the GBT Ammonia Survey, GAS; Friesen et al. 2017) find that the nearest neighbor separations of dense gas structures within the same filament are significantly larger than λ j,th and are well represented by λ cyl . The spatial resolution is limited however to approximately λ j,th ∼ 0.07 pc, and thus does not properly resolve λ j,th in sources with n 2 × 10 3 cm −3 . The above surveys support the view of hierarchical fragmentation by gravitationally unstable filaments, but lack the resolution to test what fragmentation process dominates on the scales of individual cores embedded within the clumps. The measurements of the fragmentation scale presented in §5 complement the above studies at resolutions down to ∼ 3000 au and provide further support to the view that filaments initially fragment at λ cyl and then further fragment at λ j,th .

Comparison to more active regions
Direct observations of star-forming IRDCs and embedded protoclusters have found fragmentation consistent with the thermal Jeans length (Palau et al. 2015;Beuther et al. 2015b;Teixeira et al. 2016;Busquet et al. 2016;Beuther et al. 2018) but it is unknown if these systems represent the initial state of fragmentation. Because high-mass SCCs may represent an initial stage of protocluster evolution before the formation of a high-mass star, they offer unique insight into the physical processes regulating fragmentation when compared to more evolved systems. From a survey of dense star-forming cores Palau et al. (2015) find that the fragmentation on ∼ 0.1 pc scales is best explained through thermal fragmentation. Similar results are found at sub-core spatial scales of 1000 au towards the Orion Molecular Cloud 1S (OMC-1S; Palau et al. 2017) and also consistent with the fragmentation measured in OMC-1N (Teixeira et al. 2016). While the measured median nearest neighbor separation in SCCs is consistent with the thermal Jeans length of the clump gas, the distribution also shows a distinct peak at approximately an order of magnitude higher gas density near δ nns / λ j,th ≈ 0.3 (see Figs. 12 & 13). These results may indicate continued thermal Jeans fragmentation such as in OMC-1S and OMC-1N. Beuther et al. (2015b) find results approximately consistent with thermal Jeans fragmentation towards the ∼ 800 M IRDC 18310-4, and while showing faint 70 µm emission, has similar physical properties to the SCCs in this sample. Similarly an analysis of the star-forming IRDC G14.225-0.506 favors thermal Jeans fragmentation (Busquet et al. 2016). Beuther et al. (2018) present a minimal spanning tree analysis of the separations in the CORE survey of 20 luminous (L bol > 10 4 L ) high-mass star forming regions and find fragmentation at scales on the order of the thermal Jeans length or smaller. As a possible explanation for the sub-Jeans length scales, Beuther et al. (2018) suggest that bulk motions from ongoing global collapse may have brought the fragments within closer proximity after having initially fragmented on the thermal Jeans scale. All of the sources in the CORE survey are high-mass protostellar objects (HMPOs) and more evolved than this sample. Thus our finding of fragmentation on the thermal Jeans length at an earlier evolutionary stage supports the interpretation of the COREs results and conclusion that the measured fragmentation scale may be impacted by the dynamical evolution of the protocluster.
The agreement between the nearest neighbor separations and the thermal Jeans length appears to favor a Jeans fragmentation process for stellar cluster formation. Indeed the thermal Jeans mass in typical star forming clumps is approximately 1 M , which corresponds well with the stellar mass at the peak of the Initial Mass Function. Therefore, Larson (2005) argued that the thermal Jeans process is responsible for the formation of lower mass stars in a cluster. Zhang et al. (2009) found that cores forming massive stars often 10 M , an order of magnitude greater than the thermal Jeans mass of its parental clump. These cores require additional support from turbulence to account for their formation. Furthermore, the observed measurements imply that thermal physics provide the dominant form of support, but additional models exist to describe the thermal fragmentation process that differ in geometry and density profile. For example, Myers (2017) present 2D axisymmetric models of filamentary structure that fragment through the thermal instability of Bonnor-Ebert spheres above a threshold minimum density. Because the Bonnor-Ebert radius and Jeans length have the same dependence on temperature and density with only slight differences in numerical coefficients, this leads to a fragmentation approximately equal to λ j,th . When compared to the observed median δ nns / λ j,th = 0.91 from §5, the spacings between cores predicted by Myers (2017) is δ nns / λ j,th = 0.71 (for a concentration factor q Z ≡ n /n min = 2) are broadly consistent.

Coeval formation of low-and high-mass protostars?
It is not clear if SCCs are the progenitor environments of high-mass star formation. Their high total masses (M cl ∼ 1000 M ), high central densities ( n ∼ 5 × 10 4 cm −3 ), cold gas kinetic temperatures ( T K (NH 3 ) ∼ 11 K), and low virial parameters (α vir ∼ 0.1−1) (Wienen et al. 2012;Svoboda et al. 2016) all point to persistent, bound clumps with the likely necessary physical conditions for high-mass star formation (McKee & Ostriker 2007). However, no high-mass protostars are observed. These observational facts are consistent with a scenario where high-mass stars form in SCCs through thermal fragmentation, and then accrete clump gas as initially low-mass protostars. Thus, SCCs may represent a very early and unique stage in protocluster evolution preceding the formation of high-mass protostars. This view is supported by cluster-scale theoretical simulations that incorporate protostellar and stellar feedback (Smith et al. 2009;Wang et al. 2010;Peters et al. 2010aPeters et al. ,b, 2011. Smith et al. (2009) find that no high-mass starless cores are formed in their models, and that massive stars originate from low-to-intermediate mass cores that become high-mass protostars via accretion. The mass accreted comes primarily from the surrounding clump at scales > 0.1 pc (Smith et al. 2009;Wang et al. 2010). Cyganowski et al. (2017) in a study towards the deeply embedded protocluster G11.92-0.61 discover low-mass cores in the accretion reservoir of the accreting high-mass protostellar object "MM1" with mass M * ∼ 30 − 60 M (Ilee et al. 2016). The detection of coeval low-and high-mass protostars is consistent with competitive accretion-type models of star formation (see §1). At a comparable distance of d = 3.37 +0.39 −0.32 kpc (derived from maser parallax Sato et al. 2014) and total mass to the SCCs in this study, G11.92-0.61 is more evolved, coincident with several indicators of high-mass star formation, such as Class I & II CH 3 OH masers, H 2 O masers, a GLIMPSE Extended Green Object (Cyganowski et al. 2008), numerous "hot core" molecular lines, and highvelocity collimated outflows. The sample of SCCs in this study complement the study of G11.92-0.61 in Cyganowski et al. (2017) through ALMA observations at similar resolution and sensitivity for clumps in a less active evolutionary state. In contrast, we find no clear high-mass protostellar cores or high-mass protostars in our sample of SCCs, while numerous accreting low-mass protostars are observed, as evidenced by bipolar outflows in CO/SiO. If a few of the protostars in SCCs will accrete up to high-mass stars, for which the accretion reservoir of the clump is sufficient, then these observations support a coeval mode of protocluster formation at earlier phases. When initially only low-to intermediatemass protostars are present, this coeval formation may also be termed "low-mass first" to lie in contrast to the monolithic collapse of turbulently supported high-mass cores. The competitive accretion-type simulations performed by Smith et al. (2009) find that high-mass stars form initially from intermediate mass pre-stellar cores near the center of the gravitational potential which accrete principally from collapsing clump gas up to high-mass condensations. An important feature of the Smith et al. (2009) model is that low-mass protostars form within the accretion reservoir of the central protostar, at separations < 0.15 pc. This is well matched to the distribution of nearest neighbor separations found in this work of µ 1/2 (δ nns ) = 0.118 pc (see §5). As Smith et al. (2009) point out, this signature is likely the most detectable at the early evolutionary phases of the clump where sources are less centrally concentrated in the potential and bright sources of emission are not present.
In contrast to the results of Cyganowski et al. (2017), Zhang et al. (2015) in a study of the protocluster G28.24+0.06 P1 failed to detect a distributed population of low-mass cores with Cycle 0 ALMA observations. Based on this Zhang et al. (2015) draw the conclusion that the distributed population of low-mass cores forms at a later evolutionary stage and that they are not, at least for the initial generation of protostars, coeval. Because G11.91-0.61 is at a later evolutionary stage, the distributed population of low-mass protostars observed in it may have developed after the massive cores formed. The SCCs in this study are in a similar early evolutionary phase as G28.24+0.06 and also similarly lack high-mass protostars (the maximum core mass in G28.4+0.06 is M core ∼ 16 M ). Accurate core masses are required for a quantitative analysis of the mass segregation and related length scales, but the diversity in morphologies shown within the sample, from distributed (e.g. G30660, G29558) to weakly fragmented (e.g. G28539, G29601), supports the presence of a distributed low-mass core population at the initial evolutionary phase for some systems. It is possible that depending on the initial level of support provided against fragmentation individual systems develop with varying degrees of hierarchy and segregation, and that the conclusions of Zhang et al. (2015) and Cyganowski et al. (2017) may both be correct for sources of different initial physical conditions.
The short evolutionary timescales of high-mass starless clumps, τ SCC ∼ 0.5 − 0.1 Myr for M cl = 1 − 3 × 10 3 M S16, is also consistent with the simulations of Smith et al. (2009) that show that the central, resultant high-mass protostar accretes in 0.25×t dyn ∼ 0.12 Myr the clump dynamical time, over a diameter of ∼ 0.4 pc (equivalent to the ALMA HPBW) (see also Wang et al. 2010). Similarly, (Battersby et al. 2017) perform a lifetime analysis of dense, molecular gas (N (H 2 ) 10 22 cm − 2) analyzed on a per-pixel basis from a Hi-GAL 2 deg ×2 deg field near = 30 deg. They find a timescale that is consistent for starless regions of 0.2 − 1.7 Myr, although with substantial uncertainty. The similarity in timescales is reasonable, as once a high-mass protostar forms, it would be accompanied by observational star formation indicators that identify it as a protostellar clump and remove it from the SCC category, as determined in S16. Further, we also observe hierarchical fragmentation as evidenced by the multi-modal distribution of nearest neighbor separations (see Fig. 13), as seen in G11.92-0.61. The ubiquity of filamentary structures observed (see Fig. 3) may also point to accretion mediated by sub-sonically gravitationally contracting filaments (Smith et al. 2016). This may suggest that while self-gravitating, turbulent clumps are not globally collapsing, accretion may yet be mediated through locally collapsing filaments. This latter point will be the topic of further research investigated with ALMA observations of N 2 H + J = 1 → 0 to study the kinematics of the filaments observed in this sample SCCs.

CONCLUSIONS
We present the first systematic observations of a large sample of well-vetted starless clump candidates with ALMA at high-resolution (∼ 3000 au) capable of resolving the thermal Jeans length and sensitivity (50 µJy beam −1 ) sufficient for detecting point sources down to ∼ 0.3 M and moderately compact starless cores down to ∼ 1.0 M ). The targets are selected from a complete sample of clumps identified from large Galactic Plane surveys. The sample is composed of 12 high-mass SCCs within 5 kpc from Svoboda et al. (2016) and Traficante et al. (2015) which did not show detected emission at 70 µm or other star formation indicators. Because these systems have not been affected by the extreme (proto-)stellar feedback of high-mass stars they are ideal environments to study the initial conditions of protocluster evolution. Our main findings are: 1. The newly sensitive ALMA Band 6 12 + 7 m (ν c ≈ 230 GHz) data show multiple indicators of low-/intermediate-mass star formation activity present in 11 out of 12 formerly starless clump candidates. This is determined through the presence of bipolar outflows detected in CO J = 2 → 1 and SiO J = 5 → 4 emission, and high-excitation p-H 2 CO 3 2,2 → 2 2,1 emission (E u /k = 68.1 K). These observations caution the interpretation of infrared dark clouds and SCCs identified from Galactic Plane surveys as quiescent, and unless shown otherwise are, given the findings towards this sample, likely to host low-/intermediate-mass star formation activity below the luminosity completeness of current surveys.
2. We compare representative examples of resolved and unresolved continuum sources with radiative transfer models of starless cores computed with RADMC-3D . Unresolved sources are poorly fit by starless core models with typical physical properties. The range of models does not encompass the most compact and dense cores (R flat < 1 × 10 3 au, n in 1 × 10 7 cm −3 ), but the short core free-fall times (t ff 1 × 10 4 yr) and the observed similar flux-density to Gould's Belt low-/intermediate-mass protostars, support the conclusion that these cores are protostellar even without identified outflows in CO or SiO.
3. Two high-mass starless core candidates in G28539 are identified and well fit by starless core models, with M S2 = 29 52 15 M and M S4 = 14 34 6.0 . Without supplementary measurements to infer the dust temperature profile, the masses are highly uncertain, and are consistent within the uncertainties of only forming an intermediate mass star (M * < 8 M ). 4. G28539 is the sole remaining starless clump candidate without any definitive indications of protostellar activity from the ALMA observations. It is the most massive SCC in the sample (M cl ≈ 3600 +600 −500 M , d = 4.8 +0.3 −0.3 kpc), and stands as an excellent target to study the initial conditions of protocluster evolution. A marginal 24 µm source, however, is observed coincident with 1.3 mm continuum source (G28539 S1) near the NW edge of the ALMA field, which may be evidence or protostellar activity. Further indirect evidence for star formation exists from compact SiO and CH 3 OH emission, although the source of emission is not associated with a continuum source. If these signatures are indeed associated with protostellar activity there would be no true high-mass starless clumps in this sample.

5.
A high degree of fragmentation is observed, with nearest neighbor separations consistent with the clump scale thermal Jeans length (∼ 0.1 pc). In context of previous observations that on larger scales see separations consistent with the turbulent Jeans length or cylindrical thermal Jeans length, our findings support a hierarchical fragmentation process, where the highest density regions of SCCs are not strongly supported against fragmentation by turbulence or magnetic fields.
6. Observed embedded low-to intermediate-mass star formation and thermal Jeans fragmentation in high-mass SCCs are consistent with models of star formation that form high-mass stars through gravitationally driven cloud inflow, in which low-and high-mass stars form coevally. However, further observations and followup study are necessary to properly characterize the clump star formation efficiency, protostellar accretion rates, and presence of dynamical flows in molecular tracers to validate this conclusion.

APPENDIX
A. EXAMPLE CO OUTFLOW ANALYSIS The CO J = 2 → 1 image cubes show complex emission structures that complicate the identification of coherent velocity structures such as outflows. Effects may be observed from spatial filtering, foreground and background clouds, and strong self-absorption at the clump systemic velocities. Bipolar outflows with red-and blue-shifted velocity components may still be easily observed in the data however because they are bright and are coherent in velocity over many independent channels. To illustrate these effects, we present a spatially averaged spectrum and position-velocity diagram (PV; Fig. 14) for the prominent NW-SE outflow originating from G24051 S4 (see Fig. 7). The spectrum and PV diagram are extracted from a 6. 0 diameter rectangular aperture centered along the outflow axis. Figure 14 shows bright, extended emission spanning up to ∼ 20 km s −1 from the center LSR velocity of v lsr = 83 km s −1 determined from the dense gas tracer H 2 CO 3 0,3 → 2 0,2 . The red-shifted lobe (SE) and blue-shifted lobe (NW) are clearly observed in the PV diagram at negative and positive angular offsets along the rectangular aperture axis. B. SIO J = 5 → 4 MAPS Maps of the SiO J = 5 → 4 red-and blue-shifted integrated intensities are shown in Figure 15. Three clumps have clear bipolar outflows: G24051 S5, G28565 S1, and G29601 S1. All three outflows have CO J = 2 → 1 counterparts at similar positions and velocities.
C. CORE MODEL PROPERTIES We follow a similar approach to modeling starless cores as found in Shirley et al. (2005) and Lippok et al. (2016). A similar approach is also used in McGuire et al. (2016). We assume dust opacities κ for coagulated grains and thin ice mantles in Ossenkopf & Henning (1994, hereafter OH94) for moderately processed grains with a coagulation timescale of 10 5 yr at densities between 10 4 cm −3 to 10 8 cm −3 (i.e. "OH4" through "OH6"). The coagulation density n cg from OH94 is selected for each model core based on whether the mean density (weighted by mass) is in the range 0.5 × n cg − 5 × n cg . The value of the dust opacity when interpolated at λ = 1.3 mm for 10 5 cm −3 ("OH5a") is κ = 0.90 cm 2 g −1 and varies between 0.51 cm 2 g −1 to 1.11 cm 2 g −1 over the full range of densities. We calculate total the gas mass using a dust-to-gas mass ratio of f d ≡ m d /m g = 1/110 and an ISM mean molecular weight of µ = 2.33. To fully sample the spectral range of the ISRF we extrapolate the the dust opacities from 1 µm to 90 nm using the prescription of Cardelli et al. (1989) and from 1.3 mm to 10 mm using the power law κ ν ∝ ν β with β = 1.75. In addition, scattering efficiencies for the the OH94 models are added following Young & Evans (2005) and albedos from the Weingartner & Draine (2001) WD3.1 model. The Plummer-like density profile in Eq. (4.1) is then irradiated in RADMC-3D with an external source input using the spectral energy distribution of the ISRF for a self-consistent calculation of the dust temperature distribution. We use the Black (1994) ISRF spectrum as parametrized by Hocuk et al. (2017, see Appendix B) with the UV portion of the spectrum adopted from Draine (1978). The ISRF is then varied in relative strength from the local value of the solar neighborhood by a multiplicative factor s isrf , excluding the contribution from the CMB. Figure 16a shows the ISRF specific intensity J ν for s isrf = 10 0 , 10 1 , and 10 2 with the five parametrized components clearly visible. Models are computed on a 1D radial grid from 2.5 × 10 2 au to 6.0 × 10 4 au with 100 zones with 2 × 10 6 photons to ensure convergence in the output T dust profiles over the tested range in n H . The median core mass M core integrated out a radius of 2 × 10 4 au is ∼ 1 M with the (25, 75) percentile interval ranging between 0.2 M All models Figure 16. Top: ISRF parametrization used to self-consistently calculate the temperature profiles of starless core radiative transfer models. Flux densities are scaled by factors of 1 (black), 10 1 (grey), and 10 2 (light grey), excluding the contribution from the CMB. Bottom: CDF of the gas mass enclosed within a radius of r < 2 × 10 4 au for all models (dark gray) and those with central densities nin < 3 × 10 5 cm −3 (light gray). The typical core mass is between 0.2 − 20 M .       Table 5 continued (1) Target clump name, (2) sub-structure ID number, (3) centroid right ascension coordinate, (4) centroid declination coordinate, (5) total dendrogram area, (6) Gaussian major FWHM, (7) Gaussian minor FWHM, (8) Gaussian position angle, (9) source integrated 1.3 mm flux density, (10) uncertainty in source integrated flux density, (11) source peak flux density, (12) number of bipolar CO outflows, (13) number of bipolar SiO outflows.