© ESO, 2017

1. Introduction

High-mass (O-B3 type, M_⋆> 8M_⊙) stars form in massive dense cores (MDCs, ~0.1 pc and > 10⁵ cm^-3, see Motte et al. 2007) by the accretion of gas onto stellar embryos (e.g., Beuther & Schilke 2004; Duarte-Cabral et al. 2013). The details of this process are under investigation (see recent reviews by, e.g., Tan et al. 2014; Krumholz 2015; Motte et al. 2017) and an increasing number of studies have suggested that they form through dynamical processes initiated by cloud formation (e.g., Schneider et al. 2010a; Csengeri et al. 2011; Peretto et al. 2013; Nguyen-Luong et al. 2013).

Two main theoretical scenarios have been proposed to explain the formation of high-mass stars: (1) powerful accretion driven by a high degree of turbulence (e.g., McKee & Tan 2002; Krumholz et al. 2007) or (2) colliding flows initiated by competitive accretion or cloud formation (e.g., Bonnell & Bate 2006; Hartmann et al. 2012). These two scenarios lead to distinct characteristics for the initial stages of high-mass star formation. The first model implies the existence of high-mass prestellar cores, which display high degrees of micro-turbulence. In the second model, high-mass prestellar cores never develop. When favorably located at the centers of massive reservoirs of gas, low-mass prestellar cores turn into protostars, attracting gas from further distances and eventually become high-mass protostars. Therefore, one difference between the two scenarios is the presence or absence of high-mass prestellar cores.

A decade ago, Motte et al. (2007) made the first unbiased statistical survey of MDCs. They did not find any bona-fide starless MDCs and proposed that this phase would be transitory if existent. Since then, ground-based studies have only been able to identitfy a handful of starless MDCs in molecular complexes (e.g., Russeil et al. 2010; Butler & Tan 2012) and only a few hundred starless clumps in the Milky Way (Ginsburg et al. 2012; Tackenberg et al. 2012; Csengeri et al. 2016; Svoboda et al. 2016). A few high-resolution studies have been performed with (sub)millimeter interferometers with the aim of identifying individual prestellar cores and protostars within MDCs. These attempts revealed high-mass cores that are protostellar in nature (e.g., Duarte-Cabral et al. 2013, within protostellar MDCs) or low-mass prestellar fragments (e.g., Tan et al. 2013, within starless clumps). The only prestellar core candidate of Duarte-Cabral et al. (2014, Cyg X-N53-MM2, 25 M_⊙ mass within 0.025 pc) and the single case found by Wang et al. 2014, G11P6-SMA1, ~30 M_⊙ mass within 0.02 pc) are amongst the very few good examples to date of the massive prestellar cores predicted by McKee & Tan (2002). This fact alone argues for a high-mass star formation scenario with no high-mass prestellar core phase (see Motte et al. 2017) but the debate continues.

As a first step to resolving the debate, it is crucial to find and characterize MDCs that are potential sites for the formation of high-mass stars, in molecular cloud complexes. The Herschel imaging survey of OB Young Stellar Objects (HOBYS) is dedicated to this purpose (Motte et al. 2010)¹. It is the first systematic survey of a complete sample of nearby (within 3 kpc) high-mass star progenitors at 0.1 pc scales (see initial catalogs in Nguyen Luong et al. 2011; Fallscheer et al. 2013). At five wavelengths from 70 μm to 500 μm, the HOBYS program targeted ten molecular cloud complexes forming OB-type stars, amongst which NGC 6334-6357 is one of the most massive (~7×10⁵M_⊙, see Motte et al. 2017).

NGC 6334-6357 is a molecular cloud complex that belongs to the Sagittarius-Carina arm and lies at a 1.75 kpc distance from the Sun (Matthews et al. 2008; Russeil et al. 2013). Its central part hosts a 15 pc-long filament (Russeil et al. 2013; Schneider et al. 2015), whose highest density parts consist of a ridge according to HOBYS terminology (Hill et al. 2011). NGC 6334 is itself a very active high-mass-star-forming region, as advocated by its numerous H ii regions, maser sources, and molecular outflows (see Loughran et al. 1986; Persi & Tapia 2008; Carral et al. 2002). Willis et al. (2013) identified approximately 2000 young stellar objects with Spitzer and found a heavy concentration of protostars along the NGC 6334 ridge, suggesting a mini-starburst event (in agreement with the ridge definition, see Nguyen Luong et al. 2011). Using a 1.2 mm continuum map, Russeil et al. (2010) extracted in NGC 6334 a sample of 11 clumps, including the well-known sources I and I(N) (e.g., Gezari 1982), however, only one starless candidate was identified.

In this paper, we search the NGC 6334 complex for MDCs harboring young stellar objects at the earliest phases of high-mass-star formation, that is, high-mass young protostars and high-mass prestellar cores. From Herschel and complementary images described in Sect. 2, we extracted ridges, hubs, and dense cores (see Sect. 3). We describe the SED characterization of these cores in Sect. 4. Section 5 provides a complete sample of 0.1 pc MDCs with robust mass estimates. Finally, in Sect. 6, we discuss the existence of starless MDCs and prestellar cores, the lifetimes of high-mass star formation phases, and the favored scenario for high-mass-star formation.

2. Observations

2.1. Herschel observations, data reduction, and column density images

NGC 6334 was observed by the Herschel space observatory with the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments² as part of the HOBYS (Motte et al. 2010) key program (OBSIDs: 1342204421 and 1342204422). Data were taken in five bands: 70 μm and 160 μm for PACS at FWHM resolutions or half power beam widths (HPBWs) of 5.9′′ and 11.7′′, respectively, and 250 μm, 350 μm, and 500 μm for SPIRE at FWHM resolutions of 18.2′′, 24.9′′, and 36.3′′, respectively (see Table 1). Observations were performed in parallel mode, using both instruments simultaneously, with a scanning speed of 20′′/s. Two perpendicular scans were taken, quoted as nominal and orthogonal. The size of the observed field is 2.0°× 0.5°, which corresponds to 60 pc×15 pc at a distance of 1.75 kpc. Herschel images of NGC 6334 are thus sensitive to cloud structures with ~0.1 pc to ~30 pc typical sizes, which correspond to MDCs, clumps, and clouds following definitions by, for example, Motte et al. (2017).

We reduced Herschel data using the Herschel interactive processing environment (HIPE, Ott 2010)³ software, version 10.0.2751. Versions 7.0 onwards contain a module which significantly removes the striping effects that have been observed in SPIRE maps produced with previous HIPE versions. SPIRE nominal and orthogonal maps were separately processed and subsequently combined and reduced for de-striping, relative gains, and color correction with HIPE. PACS maps were reduced with HIPE up to Level 1 and, from there up to their final version (Level 3) using Scanamorphos v21.0 (Roussel 2013). Final PACS images are shown in Figs. C.1 and C.2.

To correct for a few saturated pixels observed in SPIRE images, individual SPIRE mappings of small fields were performed with a different gain setting. The original map with saturated pixels and the small field were then combined following the method described in Appendix B of Nguyen-Luong et al. (2013) to produce full images without saturation (see Figs. C.3–C.5). Calibration for all Herschel maps was done using the offsets obtained following the procedure described in Bernard et al. (2010). This method makes use of IRIS (IRAS survey update) and Planck HFI DR2 data (Planck Collaboration I 2011) to predict an expected flux, which is compared to the median values from PACS and SPIRE maps, smoothed to the adequate resolution.

We built column density maps both at the 36.3″ and 18.2″ resolutions of SPIRE 500 μm and 250 μm data (see Figs. 1 and C.6). The procedure used to construct the 36.3″ resolution image uses the SED fit method fully described in Hill et al. (2011, 2012). The one used to construct the high-resolution column-density map is based on a multi-scale decomposition of the imaging data described in detail in Palmeirim et al. (2013). We used a dust opacity law similar to that of Hildebrand (1983) but with β = 2 instead of β = 1 and assumed a gas-to-dust ratio of 100: κ₀ = 0.1 × (ν/ 1000 GHz)² cm² g^-1.

Figure 1 presents a combined image of Herschel 70 μm and column density images of the NGC 6334 molecular complex. While the column density is concentrated along a filamentary structure, described in Russeil et al. (2013) for example, the 70 μm emission traces warm dust toward star-forming sites such as H ii regions and young stars.

2.2. Ancillary data

We complemented our Herschel observations with submillimeter and mid-infrared data (see Table 1). Multiple SCUBA-2/JCMT 450 μm and 850 μm images (Holland et al. 2013) toward the I and I(N) sources of NGC 6334 were taken from the JCMT Science Archive⁴ (Proposal ID: JCMTCAL). We reduced the data following the standard pipeline using the makemap command of the SMURF software (Chapin et al. 2013), bundled with the Gaia package. We made two mosaics from all calibrated 450 μm and 850 μm data and obtained a 0.2° × 0.2° coverage (see Figs. 1 and C.7–C.8). The angular resolutions were 8.5″ at 450 μm and 15.0″ at 850 μm. The ATLASGAL survey⁵ (Schuller et al. 2009), using the LABOCA/APEX camera at 870 μm itself covered the NGC 6334 molecular cloud with 19.2″ resolution. Dedicated SIMBA/SEST 1.2 mm observations presented by Muñoz et al. (2007) and Russeil et al. (2010) also covered the main filamentary structures of NGC 6334, with 24″ angular resolution (see coverage in Fig. 1).

At mid-infrared wavelengths, NGC 6334 was imaged by Spitzer/IRAC and MIPS at 3.6−24 μm with 1.5″−6″ resolutions, as part of the GLIMPSE and MIPSGAL surveys⁶ (Benjamin et al. 2003; Carey et al. 2009). It was also imaged by the two mid-infrared space observatories⁷ that are called WISE, at four bands and notably at 22 μm with an angular resolution of 12″, and MSX, in four bands, including 21.3 μm with spatial resolution of 18.3″ (Wright et al. 2010; Price et al. 2001).

3. Identification of dense cloud structures

We investigated the cloud structure of NGC 6334 from parsec-scale clumps (see Sect. 3.1) to 0.1 pc-scale dense cores (see Sect. 3.2).

3.1. The densest pc-scale clumps: ridges and hubs

Among parsec-scale clumps, the densest ones qualify as ridges and hubs. Hill et al. (2009) first proposed to outline ridges, that is, filamentary clumps forming high-mass stars, with a 10²³ cm^-2 column density lower limit. In contrast to threshold claims by, for example, Krumholz & McKee (2008), Hill et al. (2009) did not find any objective reason to set a definite threshold for a clump to be able to form high-mass stars. The probability distribution function (PDF) of column density in NGC 6334, studied by Schneider et al. (2015), then revealed a deviation⁸ from the first power-law tail at ~ 2 × 10²³ cm^-2, interpreted as linked to (high-mass) star formation. Figure 2 highlights the highest-density parts of NGC 6334 on a zoomed column density image. The 2 × 10²³ cm^-2 contours outline the NGC 6334 ridge spanning a 2 pc × 1 pc area and containing the I and I(N) sources as well as two smaller hubs with ~0.8 pc diameters. A narrow filament discussed at length in Matthews et al. (2008) and André et al. (2016) connects the ridge and largest hub (see Fig. 2).

We used a 2 × 10²³ cm^-2 column density background and Eqs. (1)−(2) of Nguyen-Luong et al. (2013) to derive a mass of ~15 000M_⊙ and a density of 2.5 × 10⁵ cm^-3 for the NGC 6334 ridge. With the same assumptions, the NGC 6334 hubs have masses of ~2200 and ~2800 M_⊙ and densities around 10⁵ cm^-3.

3.2. Compact sources

Compact sources were extracted using getsources (v1.131121), a multi-scale, multi-wavelength source-extraction algorithm (Men’shchikov et al. 2012; Men’shchikov 2013). Extraction is based on the decomposition of images across finely spaced spatial scales. Each of these single-scale images is cleaned of noise and background emission by iterating on the appropriate cut-off levels. Images are then re-normalized and summed over all wavelengths into a combined single-scale detection image. The main advantage of this algorithm is its multi-wavelength design. Namely, using the same combined detection image across all wavelengths eliminates the need for matching multiple catalogs obtained at different wavelengths and corresponding angular resolutions. Also, the generally non-Gaussian noise spectrum of the initial maps becomes more Gaussian on single-scale images, improving detection of the faintest objects and the deblending of sources in crowded regions. The cleaning of background emission, fully explained in Men’shchikov et al. (2012), starts from global image flattening and is iteratively refined to determine precisely the local background of each individual core.

Following the procedure established by the HOBYS consortium and fully described in this paper, the detection of compact sources is based on Herschel data plus a high-resolution column-density map. The PACS-160 μm and SPIRE-250 μm images are temperature-corrected to minimize the effects of temperature gradients. The departure to an intensity map arising from a constant temperature (17 K) is used to scale the original PACS-160 μm and SPIRE-250 μm images, with the goal of emphasizing the cold structures. The procedure is fully explained in Appendix C and images are shown in Figs. C.10−C.12. The purpose of temperature-corrected maps at 160 μm and 250 μm is to weaken the impact of external heating, especially at the periphery of photo-dissociation regions (PDRs) to help getsources focus on cold peaks rather than temperature peaks created by local heating (see Fig. C.12).

During the detection step, we simultaneously applied getsources to all images to define a catalog of sources with a unique position. At the measurement stage, we use the original (not temperature-corrected) Herschel maps from 70 μm to 500 μm plus all available submillimeter maps listed in Sect. 2.2, that is, the 450 μm and 850 μm SCUBA-2 mosaics and the 870 μm LABOCA and 1200 μm SIMBA images. These measurements go beyond simple aperture photometry since they are done together with background subtraction and deblending of overlapping sources. The resulting source catalog contains a monochromatic significance index (Sig_mono), peak and integrated fluxes (with errors), full width at half maximum (FWHM) major and minor sizes, and the position angle of the elliptical footprint for each extracted source in each far-infrared to submillimeter band. A subset of these catalogs is given in Tables A.2–A.12.

The Herschel-HOBYS imaging encompassing NGC 6334 extends from 0° to + 1.3° in Galactic latitude. On the basis of Hα, radio continuum, and HI surveys, Russeil et al. (2016) defined the velocity structure and thus the spatial area associated with the complex. Since we cannot confirm the association of NGC 6334 with gas located at the edges of the complex, we masked areas between +0.0° and +0.3° and between +1.1° and +1.3°  in Galactic latitude (see Fig. 1). In the following, we build the first-generation catalog of NGC 6334 MDCs from the corresponding getsources catalog of 4733 sources found within the unmasked area.

Sources extracted by getsources are far-infrared fluctuations among which compact dense cores persist across wavelengths. Intermediate- to high-mass dense cores⁹ are expected to be cold (10−30 K), dense cloud fragments of approximately 0.1 pc in size and several times 1−100 M_⊙ in mass. Given their expected densities, their thermal emission should mainly be optically thin for wavelengths larger than 100 μm.

Dense cores have spectral energy distributions (SEDs) peaking between 100 μm and 300 μm and should therefore have reliable flux measurements at either Herschel-160 μm or Herschel-250 μm or, more generally, at both. We define reliable flux measurements as those which have 1) signal-to-noise ratios greater than 2 for both the peak and integrated fluxes and monochromatic significance index greater than 5, and 2) deconvolved sizes less than 0.3 pc and aspect ratios less than 2. Sources in the getsources catalog which do not follow these criteria are either PDR features whose 70 μm emission is not associated with density peaks, filaments with ellipticity >2, or even ~1 pc clumps. This step selects 2063 sources out of the 4733 initial getsources sources. Roughly 90%, 6%, and 4% of the sources excluded by these criteria are 70 μm-only sources, filaments, and clumps, respectively.

160 μm is the highest-resolution Herschel wavelength (~12′′) that should correctly represent the optically thin dust emission of dense cores. In some cases, however, thermal emission at 160 μm can be contaminated by emission of small grains heated by PDRs at the peripheries of dense cores. One direct consequence of this effect is that deconvolved sizes can be larger at 160 μm than at 250 μm. The size at 160 μm or, if its flux is unreliable or contaminated, the size at 250 μm was chosen for reference to estimate the physical size of the source. Of the dense cores, 35% have their reference wavelength set to 160 μm. Another 40% of the dense cores have reliable measurements at 160 μm, but with large deconvolved sizes their reference wavelength is set to 250 μm. The remaining 25% of the dense cores only have reliable measurements at 250 μm, which thus is their reference wavelength.

A minimum number of three reliable fluxes at wavelengths larger than 100 μm are necessary to correctly constrain the SEDs of dense gas fragments. We therefore enforced a minimum of three reliable flux measurements: 1) one at either Herschel-160 μm or Herschel-250 μm, which we call the reference wavelength (see above), 2) accompanied by at least a second Herschel flux measurement at a wavelength larger than 100 μm, and 3) a third flux measurement taken at λ> 100μm with either Herschel/SPIRE, APEX/LABOCA, or SEST/SIMBA. In NGC 6334, this step removes more than half of the remaining getsources sources, leading to a sample of 654 dense cores. Sources excluded by this process are at the limit of our detection level in Herschel-SPIRE wavelengths and thus correspond to dense cores with very low mass.

4. Physical characterization of dense cores

4.1. Coherent SEDs: fluxes scaled to the same aperture

We aim to construct SEDs with fluxes measured within one aperture of the same size. Since the Herschel beam increases with wavelength, source sizes and thus the flux in an aperture often increase as well (see, e.g., Table 2 of Nguyen Luong et al. 2011). We correct for this bias by: 1) adding submillimeter fluxes with better angular resolution than our 500 μm Herschel data (e.g., 19.2′′ for APEX-870 μm and 24′′ for SIMBA-1200 μm versus 36.6′′ for Herschel-500 μm) and 2) scaling fluxes appropriately. The flux-scaling procedure primarily used by Motte et al. (2010) is described in more detail in Nguyen Luong et al. (2011), Giannini et al. (2012), and Fallscheer et al. (2013). In short, these HOBYS papers incorporated a simple linear scaling method that relies on the assumption that MDCs have intensity versus aperture radius profiles close to those found for high-mass protostellar dense cores, I ∝ θ^-1 (e.g., Beuther et al. 2002). Scaling thus provides flux estimates as they should be when integrated within the smaller aperture of the reference wavelength, that is, 160 μm or 250 μm. It aims at correcting integrated fluxes which mainly arise from optically thin thermal emission, in our case at λ> 100μm. The procedure uses deconvolved sizes, FWHM^dec, integrated fluxes, F_{λ> 100 μm}, and the formula: ${\mathit{F}}_{\mathit{\lambda }\mathit{>}\mathrm{100}\mathit{\mu }\mathrm{m}}^{\mathrm{corrected}}\mathrm{=}{\mathit{F}}_{\mathit{\lambda }\mathit{>}\mathrm{100}\mathit{\mu }\mathrm{m}}\mathrm{\times }\frac{{\mathit{FWHM}}_{\mathrm{Ref}\mathit{\lambda }}^{\mathrm{dec}}}{{\mathit{FWHM}}_{\mathit{\lambda }\mathit{>}\mathrm{100}\mathit{\mu }\mathrm{m}}^{\mathrm{dec}}}\mathrm{·}$(1)The FWHM sizes estimated using getsources are equivalent to a Gaussian FWHM, and have been deconvolved assuming that the Herschel beams themselves are Gaussian. Simple Gaussian deconvolution was stopped when sources were smaller than 0.5 × HPBW_λ, which translates into an observed FWHM$\begin{array}{c}\mathrm{obs}\\ \mathit{\lambda }\end{array}\mathrm{\le }\mathrm{1.12}\hspace{0.17em}\mathrm{\times }$ HPBW_λ. Unresolved sources were assigned at their reference wavelength a minimum physical size of FWHM$\begin{array}{c}\mathrm{dec}\\ \mathrm{160}\mathit{\mu }\mathrm{m}\end{array}\mathrm{=}\mathrm{0.5}\mathrm{\times }$HPBW_{160 μm} ≃ 5.85″ or FWHM$\begin{array}{c}\mathrm{dec}\\ \mathrm{250}\mathit{\mu }\mathrm{m}\end{array}\mathrm{\simeq }\mathrm{9.1}″$, corresponding to 0.05 pc or 0.08 pc at 1.75 kpc. Dense cores are generally resolved at their reference wavelength but 36% are not.

When fluxes remain uncorrected, SED fits of Herschel data alone lead to, on average, twice larger masses (see, e.g., Nguyen Luong et al. 2011). The reason behind this increase is that fluxes are measured within sizes varying from the reference wavelength size to approximately twice its value when measured at 500 μm. Since NGC 6334 MDCs are imaged in the submillimeter regime at higher-angular resolutions than SPIRE, adding these fluxes limits this mass increase to ~50%. As mentioned by Nguyen Luong et al. (2011), the fitted reduced χ² value is largely decreased when Herschel fluxes are linearly corrected, showing that fluxes measured in single apertures are better fitted by modified blackbody models. The linear correction is a first-order correction method of dense core fluxes measured within varying apertures. Indeed, in the case of starless dense cores, which are less centrally concentrated than protostellar dense cores, a stronger flux-scaling correction should be applied. Since there are no theoretical laws to describe the low gas concentration of a starless structure whose fragmentation is largely governed by turbulence, for consistency we decided to keep the same flux correction for all types of dense cores. In the following, the number of starless dense cores with a given mass is probably overestimated because a mass concentration steeper than Mass( <r) ∝ r (see Fig. 7) would lead to much smaller corrected 350 μm and 500 μm fluxes, higher fitted SED dust temperatures and, finally, lower derived masses.

4.2. SED fitting procedure

We used the MPFIT least-squares minimisation routine in IDL (Markwardt 2009) to fit the SEDs with modified blackbody models such as ${\mathit{S}}_{\mathit{\nu }}\mathrm{=}\mathit{A}\hspace{0.17em}{\mathit{\nu }}^{\mathit{\beta }}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathit{,}$(2)where S_ν is the observed flux distribution, A is a scaling factor, B_ν(T_dust) is the Planck function for the dust temperature T_dust, and β is the dust emissivity spectral index here set to 2.

We assumed, for the representative error of the reliable fluxes, the quadratic sum of getsources flux errors, image calibration errors¹⁰, and background extraction errors. The two latter factors combined are estimated to be of the order of 10% for 70 μm and 20% for fluxes at λ> 100 μm.

The 70 μm flux was only used for SED fitting when the temperature fitted to λ> 100 μm fluxes rises above 32 K, corresponding to a modified blackbody peaking at 100 μm. When cloud fragments are hotter than 32 K, their SEDs are incorrectly constrained by λ> 160 μm data alone. The 70 μm emission from the inner part of the protostellar envelope becomes larger than the extinction depth at 70 μm of the modified blackbody fitted at λ> 100 μm. In such cases, the 70 μm flux was therefore used as an upper limit. Throughout our sample, 20% of the sources needed the 70 μm flux in their SEDs, corresponding to sources that are either IR-bright protostellar dense cores or weak sources with poorly constrained λ> 100 μm SED portions.

Since dense cores do not have flat spectra across far-infrared to submillimeter bands, we applied color-correction factors to all Herschel fluxes. SPIRE color-correction factors depend only on the dust emissivity spectral index β, here fixed to 2. They are 4%, 4%, and 5% at 250 μm, 350 μm, and 500 μm respectively (cf. SPIRE observers manual). PACS color-correction factors themselves vary with T_dust, especially at 70 μm (Poglitsch et al. 2010). In our sample of 654 dense cores, the mean values of these corrections are 13% and 3% but they can go up to 60% and 10%, at 70 μm and 160 μm, respectively.

Since we attempted to fit all the SEDs of the 654 dense cores, even when they were very low-mass and have uncertain fluxes, we enforced fits with temperatures between 8 K and 70 K. We used the fitted reduced χ² value to discriminate between coherent sources with robust fits and inhomogeneous cloud structures, and qualified those with reduced χ²< 10 as robust SEDs, and those with reduced χ²> 10 as dubious. We found that 25% of the dense core SEDs are not robust. We excluded them in the following analysis because their parameters are too uncertain. Our fitting procedure leads to a sample of 490 dense cores for which envelope mass can be properly determined. For the remainder of the study, we consider only this sample of dense cores.

We determined the total gas + dust mass of dense cores, Mass, from their dust continuum emission, S_ν, and the following equation: $\mathit{Mass}\mathrm{=}\frac{{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}{\mathit{d}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathit{\nu }}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}}\mathit{,}$(3)where d is the distance to the Sun, 1.75 kpc, and κ_ν is the dust mass opacity. κ_ν depends on the size, shape, chemical composition, and temperature of the grains as well as the gas-to-dust mass ratio (see Ossenkopf & Henning 1994; Martin et al. 2012; Deharveng et al. 2012). κ_ν thus depends on the density and temperature of the medium considered. The main assumptions used are that the source fluxes are optically thin, sources are homogeneous in density, temperature, and optical properties along the line of sight. As in Sect. 2.1, the HOBYS consortium uses the relation κ_ν = κ₀(ν/ν₀)^β with κ₀ = 0.1 cm² g^-1 assuming a gas-to-dust ratio of 100, ν₀ = 1000 GHz and β = 2 for the dust emissivity spectral index (Hill et al. 2009; Motte et al. 2010). If we were to use β = 1.5, which is the other extreme value appropriate for dense cores, temperatures and masses would differ by 5−10%. We estimate that the overall uncertainty on masses measured through this procedure for dense cores with known distance and fixed size is approximately 50% (Roy et al. 2014).

4.3. Luminosity and volume-averaged density

We define the bolometric luminosity and volume-averaged density for our sample of dense cores. The bolometric luminosity, L_bol, should be calculated by integrating the flux densities over the source SED from 2 μm to 1 mm. To get an accurate bolometric luminosity for protostellar sources which can have significant mid-infrared emission, we used data from the GLIMPSE 3.6−8 μm, MIPSGAL 24 μm, WISE 22 μm, and MSX 21 μm catalogs¹¹. These catalogs were cross-correlated with our MDCs within a radius of 6′′, roughly corresponding to the 160 μm FWHM beam. When associations are found with both the MIPSGAL and WISE catalogs, or the MSX catalog, we favored that with the higher-resolution measurement. We estimated the mid-IR noise level in NGC 6334 as the mean dispersions in several background areas, σ_{24 μm} = 5 mJy/6″-beam.

Starless dense cores that do not show mid-infrared or even 70 μm emission can theoretically have their luminosity calculated by integrating below their respective fitted modified blackbodies. Their bolometric luminosities are thus set to the maximum value of either this estimate or the one obtained by integrating below observed fluxes.

The volume-averaged density is taken as: $\mathrm{⟨}{\mathit{n}}_{{\mathit{H}}_{\mathrm{2}}}\mathrm{⟩}\mathrm{=}\frac{\mathit{Mass}}{\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\hspace{0.17em}\mathrm{(}{\mathit{FWHM}}_{\mathrm{Ref}\mathit{\lambda }}^{\mathrm{dec}}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{3}}}\mathit{,}$(4)where $\mathit{FWH}{\mathit{M}}_{\mathrm{Ref}\mathit{\lambda }}^{\mathrm{dec}}$ is measured at the reference wavelength and deconvolved as explained in Sect. 4.1, μ = 2.8 is the mean molecular weight (Kauffmann et al. 2008), and m_H is the hydrogen mass. For unresolved sources at 160 μm or 250 μm, their upper limit size of 5.85′′ (resp. 9.1′′) leads to lower limits for their volume-averaged density defined by: $\mathrm{⟨}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{\_}\min }\mathrm{⟩}\mathrm{=}\frac{\mathit{Mass}}{\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\hspace{0.17em}\mathrm{(}\mathrm{0.5}\hspace{0.17em}\mathrm{\times }{\mathit{HPBW}}_{\mathrm{Ref}\mathit{\lambda }}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{3}}}\mathrm{·}$(5)Table 2 summarizes the main properties of our sample of 490 dense cores and Fig. 3 (top) displays their mass distribution.

5. Complete sample of NGC 6334 MDCs

The masses of dense cores extracted from Herschel data range from 0.3M_⊙ to approximately 1000 M_⊙. If we expect high-mass stars to form in the most massive dense cores, we must establish the lowest mass required by a dense core to form at least one high-mass star.

5.1. MDCs with known high-mass star formation signposts

For this purpose, we used two signposts of on-going high-mass star formation: Class II CH₃OH maser and compact radio centimeter sources. Masers result from the interaction of protostars with their envelope. In particular, 6.7 GHz Class II CH₃OH masers are found to be exclusively associated with high-mass star formation (Breen et al. 2013). As for compact radio centimeter sources, they generally arise from ionizing UCH ii regions. Dense cores displaying one of these signposts have already started to form high-mass stars. To pinpoint these high-mass star-forming dense cores, we cross-correlated the positions of our sample of 490 dense cores with those of 6.7 GHz Class II CH₃OH masers having appropriate velocities and compact radio centimeter sources listed in published catalogs (e.g., White et al. 2005; Combi et al. 2011; Pestalozzi et al. 2005; Caswell et al. 2010).

Allowing up to 10′′ offsets, we found six associations. MDC #6 is associated with both 6.7 GHz CH₃OH maser emission (Caswell et al. 2010) and an UCH ii region observed at 20 cm (White et al. 2005), while MDCs #1 (source I), #10, #14 (source I(N)), and #46 are only associated with 6.7 GHz CH₃OH maser sources (Pestalozzi et al. 2005; Caswell et al. 2010). In addition, we used N₂H⁺ observations published in Foster et al. (2011) and Russeil et al. (2010), and checked that dense cores indeed contain dense gas and that their velocities are close to the mean velocity measured for the NGC 6334 cloud complex of − 4 km s^-1, with a velocity gradient (Kraemer & Jackson 1999) from the center (V_LSR ~ 0 km s^-1) to the edge of the central filament (V_LSR ~ −11 km s^-1).

5.2. Minimum mass for MDCs in NGC 6334

To determine the minimum mass required for dense cores to form high-mass stars, we followed the method of Motte et al. (2007) and compared the dense cores of Table 2 with and without signposts of high-mass star formation.

Figure 4 shows the 24 μm flux as a function of the mass of each source, unless 24 μm flux is saturated, in which case the 22 μm or 21 μm flux is shown instead. The dense cores associated with CH₃OH masers or compact radio centimeter sources are well separated from the others. Since those dense cores with on-going high-mass star formation exhibit 75−1000 M_⊙ masses, we fixed 75 M_⊙ as the mass threshold for dense cores able to form high-mass stars in NGC 6334.

If one wants to compare the NGC 6334 mass threshold to that obtained in other studies, one must realize that the threshold depends on the typical size of the sources considered. For Cygnus X MDCs that have 30% smaller mean sizes than NGC 6334 MDCs, Motte et al. (2007) found a mass limit of 40 M_⊙ (see their Fig. 7), 45% lower than the one here. In contrast, Csengeri et al. (2014) used a much larger threshold, >650 M_⊙, to select, from the ATLASGAL survey, massive clumps with sizes of 0.4−1 pc, almost ten times larger than MDCs. The threshold found for NGC 6334 MDCs is also consistent with values used to distinguish intermediate-mass dense cores in three other HOBYS regions (Motte et al. 2010; Nguyen Luong et al. 2011; Fallscheer et al. 2013).

5.3. MDC characteristics

Using our lower mass limit of 75 M_⊙ and keeping robust dense cores with reduced χ²< 10, we ended up with a selection of 46 MDCs in NGC 6334. Table 3 lists each of these 46 MDCs along with their main physical properties as derived from the SED analysis of Sect. 4 and presented in Figs. B.1−B.46. As given in Table 4, MDCs have median FWHM sizes of ~0.08 pc, dust temperatures of ~17 K, total masses around 120 M_⊙, and median densities of 6 × 10⁶ cm^-3. Figure 3 displays their size, temperature, mass, and volume-averaged density distributions, with their mass distribution also shown on top of that of the complete sample of 490 dense cores. The role of MDCs in the high-mass star formation scenario is to be the spatial and kinematical link between ridges or hubs and individual protostars.

For the most extended MDCs, we checked for infrared associations out to a radius of 11″ from the peak. When a mid-infrared counterpart was visible but had not been measured by earlier infrared surveys, we performed aperture photometry directly on the images using the GLIMPSE and MIPSGAL guidelines for point sources¹². The resulting bolometric luminosities of NGC 6334 MDCs range from 10 L_⊙ to 9 × 10⁴L_⊙ (see Table 4).

We checked that there is a good agreement between the location of Herschel MDCs of Table 3 and the massive clumps extracted at 1.2 mm by Russeil et al. (2010). The MDCs found here are four times smaller in size and four times less massive than those clumps, that is, ~0.1 pc and 150 M_⊙ versus ~0.4 pc and 600 M_⊙, respectively. Most (8 out of 10) of the 1.2 mm clumps are subdivided into several MDCs in our data, showing that the present study allows us to follow the cloud-fragmentation process to smaller scales. The two 1.2 mm clumps which do not contain any Herschel MDCs are actually associated with compact H ii regions. Hence, their dense cloud structures could already be ionized. The masses of these two clumps were also overestimated by Russeil et al. (2010) since they underestimated their dust temperature and free-free contribution to the 1.2 mm emission. MDCs of Table 3 are generally located within the ≥ 200 M_⊙ clumps of Russeil et al. (2010), with the exception of a few that are found either close (<0.3 pc) to them or within intermediate-mass (100−200 M_⊙) clumps of Russeil et al. (2010). Therefore, the present sample of Herschel MDCs in NGC 6334 is consistent with the previous study done at a lower angular resolution and a single wavelength by Russeil et al. (2010).

6. Discussion

We now discuss the nature of the NGC 6334 MDCs (see Sect. 6.1) and use this sample to investigate the probability of the existence of high-mass prestellar cores (see Sect. 6.2), and give constraints on the high-mass star formation process (see Sects. 6.3, 6.4).

6.1. Nature and evolutionary status of MDCs

To define the nature of MDCs, we use their properties as derived from our extraction and SED fits as well as their appearance on Herschel and other submillimeter images (see Table 3 and Figs. B.1−B.46). Table 5 lists the 46 MDCs along with their SED properties and associated high-mass star formation signposts, which allow for determination of their probable nature. Figure 5 locates MDCs of all types (starless, IR-quiet, IR-bright, undefined) on the column density image of the NGC 6334 molecular complex.

6.1.1. Undefined cloud structures, doubtful MDCs

In principle, MDCs are local gravitational potential wells and should correspond to centrally concentrated cloud fragments, and thus far-infrared/(sub)millimeter emission peaks. getsources, the compact source extraction software we used (see Sect. 3.2), is powerful in deblending clusters of sources. It tends, however, to identify secondary sources around dominant ones. While some of these likely represent local gravitational wells, others may just be unbound cloud structures like the gas flows toward hubs/ridges. We therefore inspected, by eye, the column density, 160 μm, and ground-based (sub)millimeter images (see Figs. B.1−B.46) that have better angular resolution than the Herschel 250−500 μm bands. From this inspection, we confirm 32 MDC candidates and identify 14 as not centrally peaked, undefined cloud structures. Submillimeter observations at even higher-angular resolutions toward some undefined cloud structures in particular confirm that they are indeed not centrally concentrated (André et al., in prep.).

The undefined cloud structures often have SED fits of poor quality with reduced χ²> 5, suggesting they are rather diffuse. Furthermore, MDCs identified as undefined cloud structures are mostly located close to major column density peaks (see Fig. 5), whose intense far-infrared emission could have affected the source extraction. Despite our reservations, however, in the absence of a definite argument to remove these 14 undefined cloud structures from our sample, we keep them in Tables 5 and 3. Since their nature and evolutionary status are difficult to define, however, we exclude these undefined cloud structures in our following discussion.

6.1.2. Separation between IR-quiet and IR-bright MDCs

To determine the evolutionary statuses of the 32 remaining MDCs, we used the definition of IR-bright versus IR-quiet protostellar MDCs of Motte et al. (2007), followed by Russeil et al. (2010) and Nguyen Luong et al. (2011). Based on the predicted mid-infrared emission of a B3-type stellar embryo, the minimum fluxes of IR-bright MDCs located at 1.75 kpc should be S_{21 μm} = 10 Jy (Motte et al. 2010), S_{22 μm} = 12 Jy, and S_{24 μm} = 15 Jy (Russeil et al. 2010). With these criteria, IR-bright MDCs should already host at least one stellar embryo of mass larger than 8 M_⊙ while IR-quiet MDCs only contain stellar embryos of lower mass, which could further accrete. IR-quiet MDCs could therefore evolve into IR-bright MDCs of similar masses. For context, the HMPOs of Beuther et al. (2002) are IR-bright cloud structures ranging from 1 pc clumps to 0.01 pc cores, among which are several tens of IR-bright MDCs. In Fig. 4, we separate IR-bright and IR-quiet MDCs using the mid-infrared flux limits given above. The Spitzer 24 μm fluxes were primarily used but, in saturated areas, the 22 μm WISE or 21 μm MSX fluxes were the alternative choices (see Table A.8). Among MDCs of Table 5, six qualify as IR-bright protostellar MDCs: MDCs #1, #6, #10, #32, #38, and #46, as seen in Fig. 4. Four of these are also associated with CH₃OH maser sources (Pestalozzi et al. 2005).

Given their spatial sizes (~0.1 pc), MDCs cannot be precursors of single high-mass stars. MDCs with weak or even undetected mid-infrared emission probably host young protostars and/or prestellar cores. MDCs are qualified as “IR-quiet protostellar” if they are associated with compact emission at 70 μm. Most of these, however, are also detected at mid-infrared (21 μm to 24 μm) wavelengths. Similarly, a detection at both 21−24 μm and 70 μm mitigates the misidentification of externally heated diffuse cores without internal stellar embryos as being IR-quiet MDCs. “Starless” MDC candidates are generally sources for which no compact emission at 70 μm emission can be distinguished above the noise level of Herschel images. In a few cases, however, they have weak 70 μm emission at the level of or below that of the SED fits to longer far-infrared and submillimeter wavelengths (see, e.g., Figs. B.17−B.32, also Motte et al. 2010).

For each MDC, we compared the ratio of the submillimeter luminosity to the bolometric luminosity, L_bol, for the IR-quiet versus IR-bright subsamples. The submillimeter luminosity, L_submm = L_{≥ 350 μm}, is defined as the integrated luminosity below the resulting SED fit between 350 μm (857 GHz) and 1200 μm (250 GHz). The submillimeter to bolometric luminosity ratio has previously been used in the low-mass regime to distinguish between Class 0 and Class I protostars. André et al. (2000) found a separation in this ratio amongst their individual protostars (size ~ 0.01–0.1 pc) with a luminosity ratio less than 1% for Class I protostars. We cannot, however, directly use the same threshold values taken to separate Class 0 and Class I protostars. Indeed, the physical sizes, stellar luminosities, and cloud environments of MDCs and low-mass protostars are different enough to justify a change of this ratio. In the intermediate- to high-mass regime, Hennemann et al. (2010) and Nguyen Luong et al. (2011) found a clear separation between IR-bright and IR-quiet MDCs in terms of L_{≥ 350 μm}/L_bol ratio. Namely, IR-bright MDCs or UCH II regions have small (<1%) ratios and IR-quiet or starless MDCs can have ratios ranging from 1% to 20%.

Figure 6 displays the submillimeter to bolometric luminosity ratio of all NGC 6334 dense cores as a function of their mass. In this diagram, IR-bright MDCs, as defined above, are well separated from IR-quiet MDCs and lower-mass dense cores. They exhibit high masses (M ≥ 75 M_⊙) and low submillimeter to bolometric luminosity ratios (L_submm/L_bol< 0.5−1%). MDCs with L_submm/L_bol ratios larger than 1% with β = 2 or 1.5, are all IR-quiet. Among the five MDCs with L_submm/L_bol = 0.5−1%, the one associated with compact infrared emission at 8 μm and 22 μm, MDC #46, is an IR-bright protostellar MDC while the others are IR-quiet protostellar MDCs. NGC 6334-I(N) or MDC #14 itself has a L_submm/L_bol = 0.5% ratio but its mid-infrared emission is small enough for it to be an IR-quiet protostellar MDC, in agreement with the Sandell (2000) and Brogan et al. (2009) classifications of this object. Therefore, we confirm that there is a good but not perfect agreement between evolutionary classifications made with mid-infrared fluxes and L_submm/L_bol ratios.

The present sample of MDCs can help us evaluate the probability that high-mass prestellar cores exist in NGC 6334. Below, we first investigate if the starless MDCs of NGC 6334 could host high-mass prestellar cores (see Sect. 6.2) and then discuss the lifetime of any prestellar phase (see Sect. 6.3).

6.2. Quest for high-mass prestellar cores in MDCs

Our Herschel study has identified a new population of MDCs; 16 starless candidates (see Table 5). The present MDC sample can help us evaluate the probability that high-mass prestellar cores do exist in NGC 6334, especially within its starless MDCs. Thanks to our careful SED analysis using multi-wavelength Herschel images, the number of starless MDC candidates is much larger than found using the SIMBA 1 mm image only by Russeil et al. (2010), that is, 16 to 1. The dust temperatures of starless MDCs have been measured to be ~15 K (from 9.5 K to 21.1 K, see Fig. 3 middle-right), almost five degrees lower than the 20 K uniform temperature assumed by Russeil et al. (2010). With 20 K instead of 15 K, mass estimates from a flux measurement at 1 mm are underestimated by ~30%. This difference alone can explain why the number of starless MDC candidates more massive than 75 M_⊙ was largely underestimated before Herschel. Indeed, if we were to suppose a uniform temperature of 20 K for each MDC, only one starless candidate, MDC #26, would be considered massive enough. Figure 7 locates the starless MDC candidates and protostellar MDCs of NGC 6334 on a mass versus size diagram, where the gas mass concentration of other protostellar and starless samples is displayed (Bontemps et al. 2010; Butler & Tan 2012; Tan et al. 2013; Duarte-Cabral et al. 2014; Wang et al. 2014). An empirical relation close to Mass ( <r) ∝ r² is found to link starless clumps/MDCs and their prestellar core candidates. With a flux-scaling correction consistent with this relation, the number of starless MDCs in NGC 6334 reduces to 13, showing that 16 is definitively an upper limit for the number of starless MDCs in NGC 6334.

That said, the natures of the 16 starless MDC candidates of Table 5 remain in question. Nine are characterized by poorly defined SEDs (with 5 <χ²< 10, see Col. 7 of Table 3), casting doubts on the quality of their flux extraction and mass estimates, and thus their nature. Moreover, all MDCs with poor fits (except MDC #21) are located against concentrated filamentary structures (see, e.g., Figs. B.24−B.32), explaining why fluxes extracted from the highest-resolution images deviate from those obtained from their SED fits, largely set by Herschel fluxes at moderate-angular resolution. Such complex structures cast doubt on the ability of at least half of the ~0.1 pc starless MDC candidates to concentrate enough mass into a ~0.02 pc high-mass prestellar core. Of note, the location of MDC #17, close to a dominant column density source, suggests it could instead be a gas flow structure attracted by the dominant source rather than a centrally-concentrated starless MDC hosting a high-mass prestellar core.

Visually speaking, the best starless MDC candidates could be MDCs #21 (though it has a poor fit), #44, #35, #11, and #39 since they are the most isolated MDCs, with the least filamentary structures in their background. Their masses range from 80 M_⊙ to 170 M_⊙ but they have among the largest deconvolved sizes of our sample, 0.2–0.3 pc (see Col. 2 of Table 3) like in the Rosette molecular cloud (Motte et al. 2010). Their extracted characteristics lead to densities of 1−5 × 10⁵ cm^-3, more typical of low-mass prestellar cores than MDCs (10⁵ cm^-3 versus 10⁶ cm^-3, see Ward-Thompson et al. 1999; Motte et al. 2007). In contrast, mass and density-wise (150−210 M_⊙, n_H2>5 × 10⁶ cm^-3), MDCs #5, #13, #16, and #17 would be the best starless MDC candidates. As mentioned above however, MDCs #13 and #16 have poorly defined SEDs and MDC #17 may not be centrally concentrated enough. In the end, MDC #5 is our best candidate starless MDC with its unresolved size, ≤0.08 pc, and high mass, ~210 M_⊙. Interestingly, these characteristics make it at least five times denser and thus more favorable to host a high-mass prestellar core than the starless MDC candidate of Tan et al. (2013).

To explore even further the probability of NGC 6334 MDCs containing high-mass prestellar cores or high-mass protostars, we used the mass versus size diagram of Fig. 7. The masses and FWHMs of NGC 6334 MDCs are compared with those of two ~0.1 pc MDC samples for which the ~0.02 pc-scale content is known (see Bontemps et al. 2010; Tan et al. 2013). The gas mass concentration of the six most massive Cygnus X MDCs follows a classical Mass( <r) ∝ r relation, using masses and sizes estimated by Motte et al. (2007) and Bontemps et al. (2010). This relation suggests these MDCs have a centrally concentrated density distribution of ρ(r) ∝ r^-2 or steeper, similar to those of MDCs or hubs on ~0.1 – 1 pc scales (e.g., Beuther et al. 2002; Didelon et al. 2015) or low-mass protostellar envelopes on 0.01−0.1 pc scales (e.g., Motte & André 2001). By contrast, the starless structures found by submillimeter surveys or within infrared dark clouds are much less concentrated and tend to follow a Mass( <r) ∝ r² relation (Butler & Tan 2012; Tan et al. 2013). Note that the only known starless MDC N40 in Cygnus-X follows a similar relation (see Bontemps et al. 2010, Fig. 7). Given their low density (1−3 × 10⁵ cm^-3), these MDC fragments at 0.02 pc may better represent intermediate- to low-mass prestellar cores. Most starless MDC candidates in NGC 6334 are in fact 1–10 times denser than IRDC fragments (see Fig. 3 bottom-right). If they follow a Mass( <r) ∝ r² relation, they could on average form 0.02 pc cores, each of 1−5 M_⊙ mass. Given that eight starless MDCs already display non-centrally peaked morphologies, and that the seven remaining MDCs are not all very dense, the maximum number of high-mass prestellar cores expected in NGC 6334 would be seven. If MDC #5 is the only good starless MDC of NGC 6334, the minimum number of high-mass prestellar cores could be one. An ALMA project targeting the starless MDCs of NGC 6334 should soon reveal if they do indeed host groups of intermediate- to low-mass prestellar cores or simply a few dominant high-mass prestellar cores.

6.3. Constraints on the high-mass star formation process

At the chosen 0.1 pc size, our large sample of 490 dense cores should be complete for >15 M_⊙ cloud structures, according to completeness simulations performed for getsources extractions in similar Herschel images (e.g., Könyves et al. 2015; Marsh et al. 2016) and the dense core mass distribution shown in Fig. 3 top. Therefore, the present sample must be complete for the >75 M_⊙ MDCs hosting massive young stellar objects and allows us to start building an evolutionary scenario for the formation of high-mass stars. While starless MDCs cannot all form high-mass stars (see Sect. 6.2), IR-quiet MDCs will probably become IR-bright and form high-mass stars. Indeed, IR-quiet MDCs generally have more powerful SiO outflows and stronger mid-IR emission than a cluster of low-mass Class 0s (Motte et al. 2007; Russeil et al. 2010).

Table 6 lists the median characteristics of each MDC evolutionary phase and provides estimates of their statistical lifetime. With the relative number of IR-bright and IR-quiet protostellar MDCs as well as starless MDC candidates with respect to OB stars in the NGC 6334 molecular complex, we can infer the statistical lifetimes of each evolutionary phase. We use the work of Russeil et al. (2013), who performed a statistical census of O-B3 ionizing stars in the NGC 6334-NGC 6357 complex. They estimated a total number of ~150 high-mass stars in NGC 6334, which encompasses ~2/3 of the NGC 6334-6357 gas mass (~7 × 10⁵M_⊙). We assume the median age of the 150 O-B3 stars of NGC 6334 to be 1.5 ± 0.5 × 10⁶ yr according to the age of the main clusters identified by Kharchenko et al. (2013), Morales et al. (2013), and Getman et al. (2014). We also use the fragmentation level of Cygnus X protostellar MDCs (Bontemps et al. 2010) and assumed that IR-quiet and IR-bright protostellar MDCs in NGC 6334 should host, on average, two high-mass protostars. As for starless MDC candidates, given our discussion in Sect. 6.2, we assume that only one to seven of them could host at most one high-mass prestellar core. Accordingly, we obtain lifetimes of 1.1 × 10⁵ yr, 1.9 × 10⁵ yr, and 1−7 × 10⁴ yr, for the IR-bright, IR-quiet, and prestellar phases of high-mass stars in NGC 6334, respectively. The complete duration of the protostellar phase, 3 × 10⁵ yr, is only twice that of previous estimates¹³ (~1.5 × 10⁵ yr according to Motte et al. 2007; Russeil et al. 2010). The main difference is the fragmentation level here taken to be 2 while Motte et al. (2007) and Russeil et al. (2010) assumed that each MDC would form a single high-mass progenitor. The complete duration of the protostellar core phase corresponds to a few free-fall times of MDCs with ~10⁶ cm^-3 densities averaged over twice their FWHM and a single free-fall time of their ~10⁵ cm^-3 parental clump (Sect. 3.1). The lifetime of the prestellar core phase itself measured to be 1−7 × 10⁴ yr could either be as short as a single free-fall time or even be nonexistent.

Other constraints arise with the new population of 16 starless MDC candidates, which we found to account for 50% of the complete sample of MDCs. This proportion is much larger than in previous MDC studies (e.g. Motte et al. 2007; Russeil et al. 2010, 2%¹⁴ and 6%). As stated in Sect. 6.2, however, many starless MDCs listed in Table 3 may not be able to host high-mass prestellar cores. The short statistical lifetime of starless MDCs is discussed by Motte et al. (2007), who suggest that high-mass prestellar cores could be dynamically evolving into protostars.

Recent studies showed that the formation of ridges and massive stellar clusters are dynamical events, which are tightly linked (e.g., Nguyen-Luong et al. 2013; Louvet et al. 2016). Indeed, both ridges and hubs are observed to form through a global free-fall collapse (e.g., Schneider et al. 2010b; Hill et al. 2011; Hennemann et al. 2012), followed by intense formation of high-mass stars quoted as mini-starbursts (e.g., Motte et al. 2003; Nguyen Luong et al. 2011; Louvet et al. 2014). The shape and velocity drift of the south-western part of the ridge (i.e., the connecting filament in Fig. 2, see also Zernickel et al. 2013; André et al. 2016) recalls sub-filaments falling onto ridges, like those around DR21 (Schneider et al. 2010b). Interestingly, the lack of MDCs within the south-western sub-filament suggests it is stretched while falling onto the gravitation potential of the NGC 6334 ridge (see MDC distribution in Fig. 5). Moreover, a burst of star formation appears to have occurred recently in the central part of NGC 6334, corresponding to its ridge (Willis et al. 2013). In the context of the kinematical studies of ridges and hubs, the lifetime estimates of the present paper imply that either high-mass prestellar cores form and evolve in a single free-fall time of ridges/hubs or that their masses continue increasing well within the protostellar phase. The prestellar phase of high-mass star formation may be a starless MDC containing a few low-mass prestellar cores, that subsequently evolves into an IR-quiet protostellar MDC hosting a low-mass protostellar core.

The main physical properties (size, temperature, and mass) of the different evolutionary subsamples exhibit clear differences, as shown in Table 6. The median deconvolved size decreases by a factor of two during MDC evolution from starless candidate to IR-quiet and finally IR-bright MDC. While this increasing compactness suggests that MDCs continue contracting, gathering material as they evolve, it also reflects our improved ability to detect cold starless MDCs at longer wavelengths with lower angular resolution. We confirmed this decreasing size trend by comparing the 250 μm deconvolved sizes of all MDCs, and IR-bright protostellar MDCs. The average temperatures themselves rise from ~15 K to ~29 K, from starless or IR-quiet protostellar to IR-bright protostellar MDCs. Protostellar MDCs also seem to have masses increasing by 30% from their IR-quiet to their IR-bright phases. More surprisingly, starless MDCs are on average a factor of four smaller in density than protostellar MDCs: ~ 1.4 × 10⁵ cm^-3 versus ~ 5 × 10⁶ cm^-3 (see also Fig. 3), respectively. Either protostars grow in stellar mass at the same time as their MDC gas reservoir condenses or starless MDC candidates identified here do not have the ability to form high-mass stars in the near future.

The spatial distribution of the different evolutionary subsamples is shown on the high-resolution column density map in Fig. 5. MDCs are distributed along massive filaments and especially at their junctions, where material may accumulate (e.g., Schneider et al. 2012). IR-bright MDCs are found at the very center of the region, within the NGC 6334 ridge and hubs (see Sect. 3.1 and Fig. 2), suggesting, once again, that filaments and protostars grow in mass simultaneously.

6.4. Proposed scenario for the formation of high-mass stars

Combining the statistical constraints of the present paper on MDCs with results on ridges/hubs leads to a scenario for the formation of high-mass stars, which is sketched in Fig. 8 (see also Motte et al. 2017). While the association between the later H ii region phase and gas ionization by the nascent OB star UV field is rather well known (Churchwell 2002; Zavagno et al. 2007, see Fig. 8, the 6th scheme), the earliest phases have only recently been unveiled.

Ridges and hubs represent the first and largest-scale phase of high-mass star formation (see Fig. 8, 1st scheme). One of the main results of the Herschel-HOBYS survey has been the recognition of ridges and hubs as hyper-massive clumps preferentially forming high-mass stars (e.g., Hill et al. 2011; Nguyen Luong et al. 2011, see also present paper). Schneider et al. (2010b), Peretto et al. (2013), and Didelon et al. (2015), among others, showed that gas was continuously inflowing onto ridges or hubs, probably like the northern part of the connecting filament toward the NGC 6334 ridge (see Fig. 2).

Ridges and hubs play a major role in driving gas flows and matter concentration from parsec scales down to the 0.1 pc and finally 0.01 pc scales of MDCs and protostellar cores (see Fig. 8, 2nd−3rd and 4th−5th schemes). While ridges form rich clusters of high-mass stars (e.g., Hill et al. 2011; Hennemann et al. 2012; Nguyen-Luong et al. 2013), hubs may just form a couple of high-mass stars (e.g., Peretto et al. 2013; Didelon et al. 2015). These structures may represent the parsec-scale reservoir, from which gas is accreted onto 0.02 pc-scale cores.

In such a dynamical scenario, ridges/hubs, MDCs, protostellar envelopes, and stellar embryos would most likely form simultaneously. A high-mass prestellar phase may therefore not be necessary for high-mass star formation. The present study demonstrates a possible lack of high-mass prestellar cores, with the largest statistics and the most accurate mass estimates so far obtained.

At this stage, the most probable process to form high-mass stars thus is associated with converging flows and does not go through any high-mass prestellar core phase (see Fig. 8). Therefore during their starless phase, MDCs could only host a couple of prestellar cores of low to intermediate mass (2nd scheme). These prestellar cores should shortly become protostars, embedded within protostellar MDCs (3rd scheme). Protostars qualify as high-mass and IR-quiet if they are embedded within MDCs (>75 M_⊙ within 0.1 pc) but their stellar embryos remain low-mass, <8 M_⊙ (4th scheme). When large amounts of inflowing gas reach these low-mass protostars, they could become high-mass IR-bright protostars with a >8 M_⊙ stellar embryo (5th scheme), H ii regions (6th scheme), and finally O-B3 type stars.

7. Conclusions

To improve our knowledge of the high-mass star formation process, we performed an unbiased study of its earliest phases in the entire NGC 6334 molecular cloud complex. As part of the Herschel-HOBYS key program, we especially aimed at finding MDCs hosting the high-mass analogs of low-mass prestellar cores and Class 0 protostars. We used Herschel far-infrared space-based and submillimeter images together with mid-infrared and (sub)millimeter ground-based data (see Table 1) to obtain a complete census of 32 MDCs, which are the most likely progenitors of high-mass stars at 0.1 pc scales. Our main findings can be summarized as follows:

• 1.

  The Herschel imaging of NGC 6334 gives acomplete view of the cloud structures ranging from0.05 pc to 30 pc, from dense cores toclouds (see Figs. 1, and C.6–C.8).In the column density image, we outlined theNGC 6334 ridge and hubs, which are ~10⁵ cm^-3-density, 1 pc clumps, believed to form high-mass stars (see Fig. 2).

• 2.

  The exhaustive dataset gathered for the well-known NGC 6334 complex has been used as a template to determine the necessary analysis steps of HOBYS studies, focused on high-mass young stellar objects. We used getsources, an extraction software dedicated to multi-wavelength datasets. getsources performed a multi-resolution analysis of Herschel and complementary images, extracted 4733 compact sources, and provided their multi-wavelength fluxes.

• 3.

  We established a detailed procedure, for which including (sub)millimeter fluxes and scaling Herschel/SPIRE fluxes to 160 μm or 250 μm FWHM sizes are crucial steps. We visually inspected all mid-infrared to millimeter images independently and carefully constrained the complete SED of each compact source (see, e.g., Figs. B.1−B.46). Using criteria developed to identify the precursors of high-mass stars, we ensured extracted sources are genuine dense cloud structures. The resulting sample contains 490 well-characterized dense cores with typical FWHM sizes of ~0.1 pc (see Table 2).

• 4.

  We cross-correlated the NGC 6334 dense core sample with catalogs of compact radio centimeter and CH₃OH maser sources, which are signposts of high-mass star formation. We determined a lower mass limit of 75 M_⊙ for high-mass star formation to occur in the NGC 6334 region and listed the multi-wavelength fluxes and physical characteristics of 32 NGC 6334 MDCs and 14 undefined cloud structures (see Tables 5–3 and Tables A.2–A.12). MDCs have temperatures ranging from 9.5 K to 40 K, masses from 75 M_⊙ to 1000 M_⊙, and volume-averaged densities from 1 × 10⁵ cm^-3 to 7 × 10⁷ cm^-3 (see Table 4 and Fig. 3).

• 5.

  We used mid-infrared catalogs and images to expand the SEDs of the MDCs beyond far-infrared to millimeter fluxes. Integrating MDC fluxes from 3.6 μm to 1200 μm leads to bolometric luminosities ranging from 10 L_⊙ to 9 × 10⁴L_⊙ (see Table 4). We separated the MDCs into IR-bright (evolved) and IR-quiet (young) protostellar using the flux of their compact counterparts at 21−24 μm (see Fig. 4). This evolutionary discriminator is consistent with a submillimeter-to-bolometric luminosity ratio limit of L_submm/L_bol ~ 0.5−1% (see Fig. 6). Starless MDC candidates themselves are cold cloud fragments which are associated with at most a weak and compact 70 μm emission, which falls on the SED fits to their ≥160 μm fluxes.

• 6.

  Using Herschel and complementary (sub)millimeter data, we were able to constrain properly SEDs with cold temperatures. As a consequence, the present study allowed the discovery of 16 starless MDC candidates (see Figs. B.1−B.32). The nature of these MDCs must be further investigated. For various reasons, we doubt their ability to concentrate gas to high density and thus form high-mass stars in the near-future (see Fig. 7). The high-mass prestellar core phase thus remains elusive and may even not exist.

• 7.

  We used the respective proportion of MDCs found in the different evolutionary subsamples (6 IR-bright, 10 IR-quiet, and 16 starless) to infer their statistical lifetimes. Assumptions on their sub-fragmentation level lead to protostellar and prestellar lifetimes for high-mass star formation of 3 × 10⁵ yr and 1−7 × 10⁴ yr, respectively. We also found a clear increase in temperature and density in MDCs from the starless phase to the IR-quiet and IR-bright phases, suggesting that while heating its parental MDC gas, the protostellar embryo, and its parental MDC, grow in mass simultaneously (see Table 6). During star formation, MDCs tend to become more compact and cluster in denser parts of NGC 6334 (see Fig. 5), underlining the link between the formation of ridges/hubs and high-mass stars.

• 8.

  The present study shows that high-mass prestellar cores are still to be found. This lack of good candidates agrees with previous cloud structure analyses made by, for example, HOBYS. It favors a scenario wherein ridges/hubs, MDCs, and high-mass protostellar embryos form simultaneously, as shown in Fig. 8.

---

Acknowledgments

We are grateful to Alexander Men’shchikov for his help and discussions on getsources. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Part of this work was supported by the ANR (Agence Nationale pour la Recherche) project “PROBeS”, number ANR-08-BLAN-0241. We acknowledge financial support from “Programme National de Physique Stellaire” (PNPS) and program “Physique et Chime du Milieu Interstellaire” (PCMI) of CNRS/INSU, France. G.J.W. gratefully acknowledges the receipt of a Leverhulme Emeritus Professorial Fellowship.

References

Appendix A: HOBYS catalogs for the 46 MDCs found in NGC 6334

We present in this appendix the flux catalogs for the 46 MDCs discussed in the text.

Catalog entries are as follows:

• 1.

  Getsources source number.

• 2.

  Name = HOBYS_J prefix directly followed by a tag created from the J2000 sexagesimal coordinates.

• Flag for the fluxes reliability and upper limit (see criteria in Sect. 3.2). “0” means that the flux is considered unreliable because the source is too weak, S^peak/σ< 2 and/or S^int/σ< 2. Its error, 2σ, is used as an upper limit on SEDs when λ ≥ 160μm and unused otherwise (see Sect. 4.2). “1” means that the flux is considered reliable. It is used for SED fits when λ ≥ 160μm and when λ ≥ 70μm for T_dust> 32 K. “2” means that the flux is not considered reliable because the source is too large, FWHM> 0.3 pc, and/or too elliptical, e> 2. Its estimated flux is used as an upper limit on SEDs when λ ≥ 160μm or when λ ≥ 70μm for T_dust> 32 K and it is unused otherwise.

• 4.

  Estimate of the peak flux and its error in Jy/beam.

• 5.

  Estimate of the integrated flux and its error in Jy. Flux scaling has not been applied.

• 6.

  Major/minor (non-deconvolved) size estimate (FWHM) of the source in arcseconds.

• 7.

  Position angle of the major axis, east of north, in degrees.

Appendix B: Multi-wavelength images and spectral energy distribution

We present in this appendix the multi-wavelength images and spectral energy distributions (SEDs) for the 46 MDCs of NGC 6334, which are discussed in the main body of the text.

Appendix B.1: IR-bright MDCs

Appendix B.2: IR-quiet MDCs

Appendix B.3: Starless MDCs

Appendix B.4: Undefined cloud structures

Appendix C: Images used for MDC extraction

We extracted compact sources from Herschel images of Russeil et al. (2013) shown in Figs. C.1–C.5 and complementary images shown in Figs. C.6–C.12.

Figure C.6 presents the high-resolution column density map made using: 1) the procedure described in Appendix A of Palmeirim et al. (2013), based on a multi-scale decomposition of the imaging data; + and 2) the SED fitting Bayesian technique proposed by Hill et al. (2011).

Figures C.7 and C.8 present two mosaics toward the I and I(N) sources of NGC 6334 obtained from calibrated 450 μm and 850 μm SCUBA-2 data (see Sect. 2.2).

The PACS-160 μm and SPIRE-250 μm images are corrected from the effects of temperature gradients (see Figs. C.10−C.12). The purpose of the temperature-corrected maps at 160 μm and 250 μm is to weaken the impact of external heating, especially at the periphery of photo-dissociation regions to help getsources focus on column density peaks rather than emission peaks created by local heating.

To create the PACS-160 μm and SPIRE-250 μm temperature-corrected maps we first derived a basic temperature map for the region, assuming a modified blackbody model with a spectral index of 2. The color temperature at each pixel was computed by the following equation (Wien approximation): ${\mathit{T}}_{\mathrm{col}}\mathrm{=}\frac{\mathit{h}}{\mathit{k}}\mathrm{\times }\frac{{\mathit{\nu }}_{\mathrm{160}\mathit{\mu }\mathrm{m}}\mathrm{-}{\mathit{\nu }}_{\mathrm{250}\mathit{\mu }\mathrm{m}}}{\ln \left(\frac{{\mathit{S}}_{\mathrm{250}\mathit{\mu }\mathrm{m}}}{{\mathit{S}}_{\mathrm{160}\mathit{\mu }\mathrm{m}}}\mathrm{\times }\mathrm{(}\frac{{\mathit{\nu }}_{\mathrm{160}\mathit{\mu }\mathrm{m}}}{{\mathit{\nu }}_{\mathrm{250}\mathit{\mu }\mathrm{m}}}{\mathrm{)}}^{\mathrm{5}}\right)}\mathrm{·}$(C.1)The resulting T_col map was smoothed using a 40″ × 40″ resolution median filter in order to trace temperature variations

at larger scale. Flux weight-maps were then built, respectively at 160 μm and 250 μm, as ratios between modified blackbody fluxes assuming a constant temperature and the color temperature of Eq. (C.1) $\mathrm{Weight}\mathrm{-}\mathrm{map}\mathrm{(}\mathrm{160}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{160}\mathit{\mu }\mathrm{m}}\mathrm{\left(}{\mathit{T}}_{\mathrm{col}}\mathrm{\right)}}{{\mathit{S}}_{\mathrm{160}\mathit{\mu }\mathrm{m}}\mathrm{(}\mathrm{17}\mathrm{K}\mathrm{)}}$(C.2)and$\mathrm{Weight}\mathrm{-}\mathrm{map}\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{250}\mathit{\mu }\mathrm{m}}\mathrm{\left(}{\mathit{T}}_{\mathrm{col}}\mathrm{\right)}}{{\mathit{S}}_{\mathrm{250}\mathit{\mu }\mathrm{m}}\mathrm{(}\mathrm{17}\mathrm{K}\mathrm{)}}\mathrm{·}$(C.3)We used a fiducial temperature of 17 K, corresponding to the median temperature of Fig. C.9. Flux weight maps are ratio maps and therefore do not depend on the N_H2 distribution. We finally multiplied the original PACS-160 μm and SPIRE-250 μm images by their weight maps as follows: ${\mathit{I}}_{\mathrm{corrected}}\mathrm{(}\mathrm{160}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}{\mathit{I}}_{\mathrm{160}\mathit{\mu }\mathrm{m}}\mathrm{\times }\mathrm{Weight}\mathrm{-}\mathrm{map}\mathrm{(}\mathrm{160}\mathit{\mu }\mathrm{m}\mathrm{)}$(C.4)and${\mathit{I}}_{\mathrm{corrected}}\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}{\mathit{I}}_{\mathrm{250}\mathit{\mu }\mathrm{m}}\mathrm{\times }\mathrm{Weight}\mathrm{-}\mathrm{map}\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}\mathit{.}$(C.5)This weighting procedure artificially increases (decreases) the intensity of structures cooler (hotter) than 17 K, thus improving the contrast of cold and compact cores, respectively, to the emission of more diffuse, heated cloud structures. Both the PACS-160 μm and the SPIRE-250 μm temperature-corrected maps are used instead of the original PACS-160 μm and SPIRE-250 μm maps for the detection step in getsources. The temperature-corrected images are shown in Figs. C.10 and C.11.

All Tables

All Figures