© ESO 2018

1 Introduction

The central questions in high-mass star formation research focus on the fragmentation properties of the initial gas clumps that ultimately result in the final clusters, and the disk formation and accretion processes around the most massive young stars within these clusters. Furthermore, related processes such as the overall gas inflow, energetic molecular outflows, and the rich chemistry in these environments are still not comprehensively understood. For detailed discussions about these topics we refer to, for example, Beuther et al. (2007), Zinnecker & Yorke (2007), Tan et al. (2014), Frank et al. (2014), Reipurth et al. (2014), Li et al. (2014), Beltrán & de Wit (2016), and Motte et al. (2018).

Since high-mass star formation proceeds in a clustered mode at distances mostly of several kpc, high-spatial resolution is mandatory to resolve the different physical processes. In addition, much of the future evolution is likely to have been set during the earliest and still cold molecular phase, so observations at mm wavelengths are the path to follow. Most high-resolution investigations in the last decade targeted individual regions, but they did not address the topics of fragmentation, disk formation, and accretion in a statistical sense. A notable exception is the fragmentation study by Palau et al. (2013, 2014) who compiled a literature sample comprised largely of intermediate- rather than high-mass star-forming regions. However, fragmentation needs to be further studied in diverse samples, recovering larger spatial scales, and including regions of higher masses, in order to test how fragmentation behaves over a broad range of properties in high-mass star-forming regions.

To overcome these limitations, we conducted an IRAM Northern Extended Millimeter Array (NOEMA) large program named CORE: “Fragmentation and disk formation in high-mass star formation”. This program covered a sample of 20 high-mass star-forming regions at high angular resolution (~0.3″−0.4″ corresponding to roughly 1000 AU at a typical 3 kpc distance) in the 1.3 mm band in the continuum and spectral line emission. The main scientific questions to be addressed with this survey are: (a) what are the fragmentation properties of high-mass star-forming regions during the early evolutionary stages of cluster formation? (b) Can we identify genuine high-mass accretion disks, and if yes, what are their properties? Are rotating structures large gravitationally (un)stable toroids and/or do embedded Keplerian entities exist? Or are the latter embedded in the former? (c) How is the gas accumulated into the central cores and what are the larger-scale gas accretion flow and infall properties? Are the high-density cores mainly isolated objects or continuously fed by large-scale accretion flows and/or global gravitational collapse? (d) What are the properties of the energetic outflows and how do they relate to the underlying accretion disks? (e) What are the chemical properties of distinct substructures within high-mass star-forming regions?

Regarding cluster formation and the early fragmentation processes, it is well established that high-mass stars typically form in a clustered mode with a high degree of multiplicity (e.g., Zinnecker & Yorke 2007; Bonnell et al. 2007; Bressert et al. 2010; Peters et al. 2010; Chini et al. 2012; Peter et al. 2012; Krumholz 2014; Reipurth et al. 2014). Furthermore, the dynamical interactions between cluster members may even dominate their evolution (e.g., Gómez et al. 2005; Sana et al. 2012). High-spatial-resolution studies over the last decades have shown that most massive gas clumps do not remain single entities but fragment into multiple objects. However, the degree of fragmentation varies between regions (e.g., Zhang et al. 2009; Bontemps et al. 2010; Pillai et al. 2011; Wang et al. 2011, 2014; Rodón et al. 2012; Beuther et al. 2012; Palau et al. 2013; Csengeri et al. 2017; Cesaroni et al. 2017). The previous data indicate that high-mass monolithic condensations may be rare, but they could nevertheless exist (e.g., Bontemps et al. 2010; Csengeri et al. 2017; Sánchez-Monge et al. 2017). Considering subarcsecond resolution, most regions do indeed fragment, but exceptions do exist. For example, our recent investigations with the Plateau de Bure Interferometer (PdBI, now renamed to NOEMA) of the famous high-mass star-forming regions NGC 7538IRS1 and NGC 7538S revealed that NGC 7538S has fragmented into several subsources at ~ 0.3″ resolution whereas at the same spatial resolution the central core of NGC 7538IRS1 remains a single compact source (Beuther et al. 2012; see also Qiu et al. 2011 for more extended cores in the environment). At an even higher angular resolution of ≤ 0.2″ or spatial scales below 1000 AU, Beuther et al. (2013) found that even the innermost structure of NGC 7538IRS1 starts to fragment. This impliesthat the scales of fragmentation do vary from region to region. Other fragmentation studies do not entirely agree on the physicalprocesses responsible for driving the fragmentation. For example, the infrared dark cloud study by Wang et al. (2014) indicates that turbulence may be needed to explain the large fragment masses. Similarly, Pillai et al. (2011) argue for two young pre-protocluster regions that turbulent Jeans fragmentation can explain their data. However, other studies like those by Palau et al. (2013, 2014, 2015) favor pure gravitational fragmentation. Similar results are also indicated in a recent ALMA study toward a number of hypercompact HII regions (Klaassen et al. 2018). In addition to the thermal andturbulent gas properties, theoretical as well as observational investigations indicate the importance of the magnetic field for the fragmentation processes during (high-mass) star formation (Commerçon et al. 2011; Peters et al. 2011; Tan et al. 2013; Fontani et al. 2016). Furthermore, radiation feedback from forming protostars is also capable of reducing the fragmentation of the high-mass star-forming region (e.g., Krumholz et al. 2007).

It is important to keep in mind that fragmentation occurs on all scales, from large-scale molecular clouds down to the fragmentation of disks (e.g., Dobbs et al. 2014; André et al. 2014; Kratter & Lodato 2016). Different fragmentation processes may dominate on different spatial scales. In the continuum study presented here, we are concentrating on the fragmentation of pc-scale clumps into cores with sizes of typically several thousand AU. Smaller-scale disk fragmentation will also be addressed by the CORE program (see Sect. 3) through the spectral line analysis of high-mass accretion disk candidates (e.g., Ahmadi et al. 2018).

The previous investigations of NGC 7538IRS1 and NGC 7538S (Beuther et al. 2012, 2013; Feng et al. 2016) can be considered as a pilot study for the CORE survey presented here. With an overall sample of 20 high-mass star-forming regions (see sample selection below) observed at uniform angular resolution (~ 0.3″−0.4″) in the 1.3 mm wavelength band with NOEMA, we can investigate how (un)typical such fragmentation properties on core scales are. Fragmentation signatures to be investigated are, for example, the fragment mass, size and separation distributions, and how they relate to basic underlying physical processes.

In this paper, we present the sample selection, the general survey strategy as well as the observational characteristics. The remaining of the paper focuses on the continuum data and fragmentation properties of the sample. The other scientific aspects of this survey will be presented in separate publications (e.g., Ahmadi et al. 2018; Mottram et al. in prep.; Bosco et al. in prep.).

2 Sample

Our sample of young high-mass star-forming regions was selected to fulfill several criteria. (a) Regions of our sample have luminosities > 10⁴ L_⊙ indicating that at least an 8 M_⊙ star is forming. (b) Distances are limited to below 6 kpc to ensure high linear resolution (~1000 AU). (c) The sources are high-declination (decl. > 24°) to obtain the best possible uv-coverage (implying that they are either not at all or at most poorly accessible with the Atacama Large Millimeter Array, ALMA). Furthermore, only sources with extensive complementary high-spatial resolution observations at other wavelengths were selected to better characterize their overall properties. In this context, the sample is also part of a large e-Merlin project led by Co-I Melvin Hoare to characterize the cm continuum emission of the sample at an anticipated spatial resolution of down to 30 mas. The initial luminosity selection was based on luminosity and color–color criteria. Figure 1 presents the corresponding luminosity–color plot. We used the luminosity-color plot as a sample selection tool as the y- and x-axes act as proxies for stellar mass and evolutionary stage, respectively. By the time massive forming stars have reached 10⁴ L_⊙, the luminosity is determined primarily by the stellar mass as at this stage the accretion luminosity only contributes a small fraction of the total luminosity even at high-accretion rates (e.g., Hosokawa & Omukai 2009; Hosokawa et al. 2010; Kuiper & Yorke 2013; Klassen et al. 2016). We also expect that over time the IR colors will evolve from red to blue as the envelope material is dispersed and/or accreted (e.g., Zhang et al. 2014).

Many sample sources are covered by the Red MSX sources survey (RMS; Lumsden et al. 2013), and a few additional prominent northern hemisphere regions are included as well. Our sample excludes the few sources that fulfill these selection criteria but which already have been observed at mm wavelengths with high angular resolution (e.g., W3IRS5, NGC 7538IRS1/S; Rodón et al. 2008; Beuther et al. 2012). The resulting sample of 18 regions is complete within these described selection criteria. Because NGC 7538IRS1 and NGC 7538S were observed in an almost identical setup (only the compact D-array data were not taken), they are considered as a pilot study and their results are incorporated into the analysis of the CORE project. Table 1 presents a summary of the main source characteristics, including their local-standard-of-rest velocity v_lsr, distance D, luminosity L, mass M (see also Sect. 5.3), their 8and 21 μm fluxes, H₂O, CH₃ OH maser and cm continuum associations as well as references for the distances and luminosities. Figure 2 shows a larger-scale overview of the twenty regions with the near- to mid-infrared data shown in color and the 850 μm continuum single-dish data (Di Francesco et al. 2008) presented in contours.

Regarding the evolutionary stage of the sample, they are all luminous and massive young stellar objects (MYSOs) or otherwisenamed high-mass protostellar objects (HMPOs). Subdividing the regions a bit further, some regions show very strong (sub)mm spectral line emission indicative of hot molecular cores (AFGL2591, G75.78+0.34, CepAHW2, W3(H₂O), NGC 7538IRS1), other regions are line-poor (e.g., S87IRS1, S106, G100.3779, G084.9505, G094.6028, G138.2957, G139.9091), and the remaining sources exhibit intermediate-rich spectral line data. Furthermore, the sample covers various combinations of associated cm continuum, H₂ O and class II CH₃ OH maser emission (Table 1). Following Motte et al. (2007), we checked whether the sources belong to the so-called IR-bright or IR-quiet categories with the dividing line defined as IR-quiet when $S_{21\mu \text{m}}<10\text{\hspace{0.17em}Jy}{\left(\frac{1.7\text{\hspace{0.17em}kpc}}{D}\right)}^2$ $\left(\frac{L}{1000L_{\odot }}\right)$. In contrast to our initial expectation that all sources would classify as IR-bright, we clearly find some diversity among the sample (see Table 1). While the majority indeed qualified as IR-bright, a few sources fall in the IR-quietcategory. Perhaps slightly surprising, a few of our line-brightest sources are categorized as IR-quiet (e.g., CepA and W3(H₂O)). Therefore, the differentiation in these two categories only partly implies that the IR-quiet sources are potentially younger, but it suggests at least that these sources are still very deeply embedded into their natal cores. In this embedded stage, they are already capable of driving dynamic outflows, have high luminosities, and produce a rich chemistry.

A different evolutionary time indicator sometimes used is the luminosity-over-mass ratio L∕M of the regions (see Table 1; e.g., Sridharan et al. 2002; Molinari et al. 2008, 2016; Ma et al. 2013; Cesaroni et al. 2017; Motte et al. 2018). The CORE sample covers a relatively broad range in this parameter space between roughly 20 and 700 L_⊙/M_⊙. However, this ratio is not entirely conclusive either. For example, the region with our lowest ratio (S87IRS1 with L∕M~ 18L_⊙/M_⊙), that could be indicative of relative youth, is classified otherwise as IR-bright which seems counterintuitive at first sight. Since the various age-indicators are derived from parameters averaged over different scales, it is possible that they are averaging over subregions with varying evolutionary stages and are hence not giving an unambiguous evolutionary picture.

In summary, the CORE sample consists of regions containing HMPOs/MYSOs above 10⁴ L_⊙ from the pre-hot-core stage to typical hot-cores and also a few more evolved regions that have likely already started to disrupt their original gas core. The evolutionary stages are comparable to the sample by Palau et al. (2013, 2015) with the difference that they had a large fraction of sources below 10⁴ L_⊙ and even below 10³ L_⊙ (only four regions above 10⁴ L_⊙).

Fig. 1 Sample selection plot where the luminosity (in units of L⊙) is plotted against the MSX 21/8 μm color. Horizontal bars mark uncertainties in the color. While the blue sources fulfill our selection criteria, the red ones are below our luminosity cut of 104 L⊙. Green sources are those for which high-resolution mm data already exist and which were therefore excluded from the observations.

Fig. 2 Large-scale overview images for the whole CORE sample. The color scale shows three-color images with blue, green, and red from Spitzer 3.6, 4.5, and 8.0 μm for all sources except IRAS 23033, IRAS21078, G100, G094, and IRAS 23385 for which WISE 3.4, 4.6, and 12 μm data are presented. Furthermore, W3IRS4 uses Spitzer 3.6, 4.5 μm, and MSX 8 μm. The contours show SCUBA 850 μm continuum data (di Francesco et al. 2007; contour levels 20%, 40%, 60%, and 80% of the peak emission) for all sources except G100, G084, and G108 where these data do not exist.

3 CORE large program strategy

Following our experience with NGC 7538IRS1 and NGC 7538S (Beuther et al. 2012, 2013), we devised the CORE survey in a similar fashion. The full sample is observed in the 1.3 mm band, and a subsample of five regions will also subsequently be observed at 843 μm. Here we focus on the 1.3 mm part of the survey for the full sample. The shorter wavelength study will be presented once it is complete.

Several aspects were considered to achieve the goals of the project: (i) the most extended A-configuration of NOEMA was used for the highest possible spatial resolution. (ii) Complementary observations with more compact configurations of the interferometer recover information on larger spatial scales. Simulations showed that adding the B and D configurations provided the best compromise between spatial information and observing time. (iii) To also cover very extended spectral line emission, short spacing observations from the IRAM 30 m telescope were added. (iv) Spectrally, among other lines our survey covers CH₃ CN to trace high-density gas as might be found in accretion disks and/or toroids (e.g., Cesaroni et al. 2007) and H₂ CO which traces lower-density, larger-scale structures. Both, CH₃CN and H₂CO are also well known temperature tracers (e.g., Mangum & Wootten 1993; Zhang et al. 1998; Araya et al. 2005). Furthermore, outflow tracers such as ¹³CO and SO are included. A plethora of additional lines are also covered to investigate the chemical properties of the regions. An early example of such investigation can be found in the paper about the pilot study sources NGC 7538IRS1 and NGC 7538S by Feng et al. (2016).

With the wide-band correlator units WIDEX, a spectral range from 217.167 to 220.834 GHz was covered at a spectral resolution of 1.95 MHz, corresponding to a velocity resolution of ~2.7 km s⁻¹ at the given frequencies. Figure 3 shows an example spectrum from AFGL2591. These wide-band units are used to extract the line-free continuum, as well as to get a chemical census of the region. Furthermore, the velocity resolution is sufficient for outflow investigations. However, to study the kinematics of the central rotating structures, higher spectral resolution is required. Therefore, we positioned the eight narrow band correlator units to specific spectral locations that covered the most important lines at a spectral resolution of 0.312 MHz, corresponding to a velocity resolution of ~0.43 km s⁻¹ at the given frequencies. Table 2 shows the spectral lines covered at this high-spectral resolution. For more details about the spectral line coverage we refer the reader to the CORE paper by Ahmadi et al. (2018).

For the complementary IRAM 30 m short spacings observations, we mapped all regions with approximate map sizes of 1′ in the on-the-fly mode in the 1 mm band. Since the bandpasses at the 30 m telescope are broader and the receivers work in a double-sideband mode, the 30 m data cover a broader range of frequencies between ~213 and ~221 GHz in the lower sideband and between ~229 and ~236 GHz in the uppersideband. The line data that are covered by the NOEMA and 30 m observations can be merged and imaged together, whereas the remaining 30 m bandpass data can be used as standalone data products. Since we did not use the single-dish data for the continuum study presented here, we refer to the CORE paper by Mottram et al. (in prep.) for more details on the IRAM 30 m data.

More details about the CORE project are provided at the team web-page¹. There, we will also provide the final calibrated visibility data and imaged maps. The data release will take place in a staged fashion. The continuum data are published here, and the corresponding line data will be provided with subsequent publications.

4 Observations

The entire CORE sample (except the pilot sources NGC 7538IRS1 and NGC 7538S) was observed at 1.37 mm between summer 2014 and January 2017 in the three PdBI/NOEMA configurations A, B and D to cover as many spatial scales as possible (see Sect. 3). The baseline ranges for all tracks in terms of uv-radius are given in Table 3. The shortest baselines, typically between 15 and 20 m, correspond to theoretically largest recoverable scales of 16″− 20″. For each track, two sources were observed together in a track-sharing mode. The phase centers of each source and the respective source pairs for the track-sharing are shown in Table 1. Since each source was observed in three different configurations, at least three (half-) tracks were observed per source. Depending on the conditions, several source pairs were observed in more than three (partial) tracks in order to achieve the required sensitivity and uv-coverage. Altogether, this multi-configuration and multi-track approach resulted in excellent uv-coverage for each source, an example of which is shown in Fig. 4. Typically two phase calibrators were observed in the loops with the track-sharing pairs. For the final phase calibration, we mostly used only the stronger ones. Depending on array configuration and weather conditions, the phase noise varied between ~10 and ~50deg. Bandpass calibration was conducted with observations of strong quasars, for example, 3C84, 3C273, or 3C454.3. The resulting spectral baselines are very good, over the broad WIDEX bandpass as well as the narrow-band bandpasses (e.g., see Fig. 3). The absolute flux calibration was conducted in most cases with the source MWC349 where an absolute model flux of 1.86 Jy at 220 GHz was assumed². For only very few tracks in which that source was not observed, the flux calibration was conducted with other well-known calibrators (e.g., LKHα101). The absolute flux scale is estimated to be correct to within 20%.

To achieve the highest angular resolution, uniform weighting was applied during the imaging process. The final synthesized beams for the continuum combining all NOEMA data vary between ~ 0.32″ and ~ 0.5″ with exact values for each source given in Table 3. The full width at half maximum of the primary beam of our observations is ~22″. To create the continuum images, we carefully inspected the WIDEX bandpasses for each source individually and created the continuum from the line-free parts only. The 1σ continuum rms correspondingly varies from source to source. This depends not only on the chosen line-free channels, but also on the side-lobe noise introduced by the strongest sources in the fields. Although the uv-coverage is very good(Fig. 4), not all side-lobes can be properly subtracted, and the final noise also depends on the uv-coverage.

To reduce calibration, side-lobe and imaging issues, we explored how much self-calibration would improve the data quality. For that purpose, we exported the continuum uv-tables to CASA format and did the self-calibration within CASA (version 4.7.2; McMullin et al. 2007). We performed phase self-calibration only, and the time intervals used for the process varied from source to source depending on the source strength. Solution intervals of either 220, 100, or 45 s were used, where 45 s is the smallest possible interval due to averaging of the data during data recording. Interactive masking during the self-calibration loops was applied, with only the strong peaks used in the first iterations and then subsequently adapted to the weaker structures. After the self-calibration, we again exported the data to GILDAS format and conducted all the imaging within GILDAS to enable direct comparisonswith the original datasets. Again, uniform weighting was applied and we cleaned the data down to a 2σ threshold. To show the differences of the images prior to and after the self-calibration process, Appendix C presents thederived images before and after the self-calibration. The contouring is done in both cases in 5σ steps. Careful inspection of all data shows that no general structural changes were created during the self-calibration process. The self-calibration improved the data considerably with reduced rms noise and slightly increased peak fluxes. We find that the flux-ratios between the main substructures within individual regions remained relatively constant prior to and after self-calibration. Below, we conduct the analysis with the self-calibrated dataset. Table 3 presents the 1σ continuum rms for all sources before and after self-calibration. We typically achieve sub-mJy rms with a range between 0.05 and 1.9 mJy beam⁻¹ for the 18 new targets. Only the pilot source NGC 7538IRS1 has a slightly higher rms of 5.2 mJy beam⁻¹ which can be attributed to the higher source strength and the missing D-array observations. Primary-beam correction was applied to the final images, and the fluxes were extracted from these primary-beam corrected data (Sect. 5.2). Evaluating the measured peak flux densities S_peak and noise values (rms_sc) in Table 3, we find S/N of between 39 and 326 with the majority of region (13) exhibiting S/N greater than 100. We provide the original pre-self-calibration images, the images after applying self-calibration as well as the primary-beam corrected images in electronic form³.

Fig. 3 Example wide-band spectrum extracted toward AFGL2591. The most important lines in the bandpass are marked.

Fig. 4 Example uv-coverage for CepA. The different colors correspond to different observed (half-)tracks. Red and black correspond to D-array observations, blue and cyan to B-array, and green to A-array data.

Simulated observations

To better understand how the imaging affects our results, we simulated a typical observation. The details of the simulations can be found in Appendix B. In the following, we summarize the method and results. We used real single-dish dust continuum data from the large-scale SCUBA-2 850 μm map of Orion by Lane et al. (2016), converted the flux to 1.37 mm wavelength (assuming a ν^3.5 frequency-relation), rescaled the spatial resolution and flux density to a distance of 3 kpc, and imaged different parts of Orion with the typical uv-coverage and integration time from the CORE project. Similar to our observations, the rms varied depending on whether a strong source (in this case Orion-KL) was present in the observed field. While the point source mass sensitivity is very good, between 0.01 and 0.1 M_⊙ (depending on the rms), with our spatial resolution typical Orion cores are extended structures, rather than point sources, even at a distance of 3 kpc. Hence, the dependence of the rms noise on the strongest sources in the field strongly affects the actual core mass sensitivity for extended structures as well. Taking the two examples shown in Appendix B, cores with masses down to ~1 M_⊙ are detectable in fields without very strong sources. If such a low-mass core were within the stronger Orion-KL field, it would not be detectable anymore. Therefore, the core mass sensitivities strongly depend on the strongest and most massive sources within the respective observed fields. The dynamic range limit of the simulations of Orion-KL is approximately 53.

5 Continuum structure and fragmentation results

5.1 Source structures

Figures 5 and 6 present the 1.37 mm continuum data of the full CORE sample. While Fig. 5 shows the data in angular resolution over the full area of the primary beam of the observations, Fig. 6 uses the distances of the sources (Table 1) and presents the data at the same linear scales, making direct comparisons between sources possible. The first impression one gets from these dust continuum images is that the structures are far from uniform. While some sources are dominated by single cores (e.g., IRAS 23151, AFGL2591, S106, NGC 7538IRS1), other regions clearly contain multiple cores with a lot of substructures (e.g., S87IRS1, IRAS 21078, W3IRS4), some of which have more than ten cores within a single observed field (see Sect. 5.2). We see no correlation between the number of fragments and the distances to the sources. We discuss this fragmentation diversity in more detail in Sect. 6.

Fig. 5 Compilation of 1.37 mm continuum images for CORE sample on the same angular scale. The contouring is in 5σ steps (see Table 3). The sources are labeled in each panel, and the synthesized beams are shown at the bottom-left of each panel. A comparison figure converted to linear scales is shown in Fig. 6. Magnifications and absolute flux-scales are shown in Appendix C.

Fig. 6 Compilation of 1.37 mm continuum images for CORE sample converted to linear resolution elements. The contouring isin 5σ steps (see Table 3). The sources are labeled in each panel, and the synthesized beams are shown at the bottom-left of each panel.

5.2 Source extraction

To extract the sources from our 20 images, we used the classical CLUMPFIND algorithm by Williams et al. (1994) on our self-calibrated images. As input parameters we used the 5σ contour levels presented in Figs. 5 and 6 as well as in Appendix C. These images sometimes also show negative 5σ contours, indicating that the interferometric noise is neither uniform nor really Gaussian. Therefore, we inspected all sources identified by CLUMPFIND individually and only included those where the peak flux density is ≥10σ (two positive contours minimum in Appendix C). The derived positional offsets from the phase center, peak flux densities S_peak, integrated flux densities S_int and equivalent core radii (calculated from the measured core area assuming a spherical distribution) are presented in Table A.1 (S_peak and S_int are derived from the primary-beam corrected data).

To estimate the amount of missing flux filtered out by the interferometric observations, we extracted the 850 μm peak flux densities from single-dish observations, mainly from the SCUBA legacy archive catalog (Di Francesco et al. 2008). Since this dataset has a final beam size of 22.9″ it covers our primary beam size very well. Scaling this 850 μm data with a typical ν^3.5 dependency to the approximate flux at our observing frequency of 220 GHz, we can compare these values to the sum of the integrated fluxes measured for each target region from our previous CLUMPFIND analysis. Table 3 presents the corresponding missing flux values (mf in percentage) for the sample (for two regions – G100 & G108 – we did not find corresponding single-dish data). The amount of missing flux varies significantly over the sample, typically ranging between 60 and 90%. The only extreme exception is W3(H₂O) where only 25% of the flux is filtered out. This implies that for this region the flux is strongly centrally concentrated without much of a more extended envelope structure. For the remaining sources, even with the comparably good uv-coverage (Fig. 4) a significant fraction of the flux is filtered out. The variations from source to source indicate that the spatial density structure also varies strongly from region to region (see discussion in Sect. 6).

There is a broad distribution in the number of cores identified in each region. We find between 1 and 20 cores among the different regions (see Table 4). To check whether this range of identified cores is related to our mass sensitivity, in Fig. 7 we plot the 5σ mass sensitivity (Table 3 and Sect. 5.3) vs. the number of identified cores (excluding NGC 7538IRS1 because of its unusually poor mass sensitivity limit, Table 3). While there might be a slight trend of more cores toward lower mass sensitivity limits, our lowest mass sensitivity limit region CepA also shows only two cores. In the main regime of 5σ mass sensitivities between 0.1 and 0.3 M_⊙ we do not see a relation between the number of identified cores and the mass sensitivity. Hence, the number of identified cores does not seem to be strongly dependent on our mass sensitivity limits below 0.4 M_⊙.

Fig. 7 Number of identified cores plotted against the 5σ mass sensitivity. NGC 7538IRS1 is excluded because of its unusually high mass sensitivity limit.

5.3 Mass and column density distributions

Assuming optically thin dust continuum emission at 220 GHz, we can estimate the gas masses and peak column densities for all identified cores in the sample. Following the original outline by Hildebrand (1983) in the form presented by Schuller et al. (2009), we use a gas-to-dust mass ratio of 150 (Draine 2011), a dust mass absorption coefficient

κ of 0.9 cm² g⁻¹ (Ossenkopf & Henning 1994 at densities of 10⁶ cm⁻³ with thin ice mantles) and average temperatures for each region derived from the CORE IRAM 30 m H₂ CO data. H₂ CO is a well-known gas thermometer in the ISM (Mangum & Wootten 1993), and we derive beam-averaged temperatures from the single-dish spectra toward the peak positions of each region at a spatial resolution of 11″. For the temperature estimates we fitted the data with the XCLASS tool (eXtended CASA Line Analysis Software Suite) tool (Möller et al. 2017). XCLASS models the spectra by solving the radiative transfer equation for an isothermal homogeneous object in local thermodynamic equilibrium (LTE), using the molecular databases VAMDC and CDMS⁴ (Müller et al. 2001). XCLASS employs the model optimizer package MAGIX (Modeling and Analysis Generic Interface for eXternal numerical codes) to find the best fit solutions (Möller et al. 2013). The derived temperatures are shown in Table 3. Since we derive beam-averaged temperatures from the single-dish data, the actual temperatures of individual cores at smaller spatial scales may vary compared to those beam-averaged temperatures. More detailed temperature analysis from the combined interferometer plus single-dish data is beyond the scope of this paper and will be conducted in future work on the CORE data. The mass estimates are in general lower limits since we are filtering out large-scale flux that may be associated with the dense cores (see also Appendix B). Furthermore, while the optically thin assumption for the dust emission should be valid in most cases, there may be some exceptions like CepA where high peak flux densities (Tables 3 and A.1) imply high brightness temperatures indicating moderate optical depth at these peak positions. However, since the masses are calculated typically over areas larger than just the peak, and the brightness temperatures decrease quickly with distance from the peak, this effect should be comparably weak.

The derived core masses and column densities are presented in Table A.1 and roughly span 0.1 to 40 M_⊙, and 5 × 10²² to 10²⁵ cm⁻². For the mass and column density analysis, we excluded sources for which the continuum emission is clearly dominated by HII regions and hence show barely dust continuum emission. These are specifically W3(OH) (cores #1 and #2 in W3(H₂O), the southern ring-like region in W3IRS4 (sources #5 and #6) and core #2 in S87IRS1. For several other cores, the fluxes were corrected for free–free emission for the mass determinations (see Table A.1).

Using similar assumptions, we also re-estimated the large-scale mass reservoir for the sample. For most sources, we used the integrated 850 μm fluxes derived by Di Francesco et al. (2008), while for IRAS 23151 the 1.2 mm flux was derived from the MAMBO data presented in Beuther et al. (2002), and for G084 we used the 1.1 mm BOLOCAM data from Ginsburg et al. (2013). The used gas-to-dust mass ratio and average H₂ CO derived temperatures are the same as above, and we used for the single-dish data dust absorption coefficients κ of 0.78, 0.9, and 1.4 cm² g⁻¹ at 1.2, 1.1, and 0.85 mm wavelengths, respectively (Ossenkopf & Henning 1994; at densities of 10⁵ cm⁻³). The derived total masses are presented in Table 1 (for G100 and G108, the masses are taken from C¹⁸ O(3–2) data from Maud et al. 2015). While the regions have typical mass reservoirs of several 100 M_⊙, the sample spans a comparably broad range between ~40 and ~1500 M_⊙ (for G108 even higher masses are measured, however over a comparably large area with radius 1.4 pc, Table 1 and Maud et al. 2015).

For the NOEMA-only derived core parameters, Fig. 8 shows histograms of the masses and column densities. The combined mass distribution shows that most detected cores are in the range between ~0.1 and ~10 M_⊙ with only a few cores exceeding 10 M_⊙. The most massive coreis in NGC 7538IRS1 with 43 M_⊙ (although significant free-free contamination may affect the estimate for this source, Beuther et al. 2012). Regarding the cores in excess of 10 M_⊙, there is no clear trend whether they are found as isolated objects or embedded in fragmented regions. For example, comparably massive cores are found in the low-fragmentation regions NGC 7538IRS1 or AFGL2591, but cores of similar mass are also found in more fragmented regions like IRAS 23151, IRAS 23033, G75.78, as well as in the intermediately fragmented region W3(H₂O). The peak column densities are very large, typically exceeding 10²³ cm⁻² and even going above 10²⁵ cm⁻² for a few exceptional regions. Figure 9 plots the column densities against the masses, and while we see a scatter, there remains nevertheless a trend that column densities and masses are correlated. If one takes into account the distance-dependencies of our derived parameters (color-coding in Fig. 9), we see that the higher-mass-lower-column-density sources are found at on average larger distances. With increasing distance, the physicalsize of the beam (where the column density is measured within), increases as well. Such larger area beams coverthe central highest-column-density peak position but also more lower-column-density environmental gas. This smoothing slightly decreases the measured column densities with increasing distance. The other way round, increasing the covered area with distance also increases the measured masses. Hence, part of the scatter in Fig. 9 is caused by the distance range of our sample. For smaller distance bins, the scatter is significantly reduced.

Using the derived equivalent radii of the cores from the CLUMPFIND analysis (Table A.1), we can also derive mean densities for all cores under the assumption of spherical symmetry. Figure 10 plots these mean densities against the corresponding core masses, again color-coded with distance. While these average densities are rather high, typically between 10⁶ and 10⁸ cm⁻³, there is no clear trend between the densities and the masses. Taking the distances into account again, the scatter is reduced but identifying trends within distance-limited ranges is still difficult. Hence, in this sample, the core densities are similar over the whole range of observed core masses. Having a correlation between mass and column density but less good correlation between mass and average density implies that the core masses should correlate with their sizes (their equivalent radii). Figure 11 presents the corresponding data again color-coded with distance. And indeed mass and size are well correlated for the sample, again much tighter if one looks at limited distance ranges. Figure 11 also plots lines of constant column densities between 10²³ and 10²⁵ cm⁻². While most regions scatter between the 10²³ and 10²⁴ cm⁻² lines, also subsamples between limited distance ranges do not follow constant column density distributions but increase in column density withincreasing mass, as already shown in Fig. 9.

To estimate the typical Jeans fragmentation lengths and masses for the clump scales, we assume mean densities of the original larger-scale parental gas clumps between 10⁵ and 10⁶ cm⁻³ (e.g., Beuther et al. 2002; Palau et al. 2014) and a temperaturerange between 20 and 50 K, typical for regions in the given evolutionary stages. For such conditions, the estimated Jeans length is between ~5500 and 27 700 AU. For comparison, the corresponding Jeans masses in this parameter range vary between 0.3 and 3.5 M_⊙. While a large fraction of the core masses lies within the regime of the Jeans masses, a non-negligible number of sources also have higher masses (~36%) in excess of the Jeans mass of the original cloud. Since our mass estimates are lower limits, in reality even more cores may exceed the estimated Jeans masses. However, since the mass estimates are affected by many uncertainties (in addition to the missing flux, the assumeddust properties and temperatures are adding an uncertainty of factors between two and four), the core separations may be a better proxy for analyzing the fragmentation properties of the gas clumps.

Fig. 8 Histograms of masses (top panel) and column densities (bottom panel) for all detected cores.

Fig. 9 Gas column densities vs. masses for all detected cores. The color-coding shows the distances of the sources.

Fig. 10 Mean core densities vs. masses for all detected cores. The color-coding shows the distances of the sources.

Fig. 11 Core masses vs. equivalent radii for all detected cores. The color-coding shows the distances of the sources. The dashed lines show constant column density with levels of 1023, 1024, and 1025 cm−2 from right to left.

5.4 Core separations

To quantify the core separations in all 20 sample regions, we employed the minimum spanning tree algorithms available within the astroML software package (VanderPlas et al. 2012) which determines the shortest distances that can possibly connect each of the cores in the sampled field. From this, the minimum, maximum, and mean separations of the cores in each field were determined, and are presented in Table 4, with the distribution of nearest neighbor separations shown in Fig. 12. Since our data are 2D projections of 3D distributions, these measured separations are necessarily lower limits. The minimum core separations are typically on the order of a few 1000 AU (peak at ~2000 AU, similar to Palau et al. 2013) with only a few core separations for the most nearby sources being measured below 1000 AU. However, this lower limit is most likely not a real physical lower separation limit but associated with the spatial resolution. With typical resolution elements around 0.3″−0.4″ (Table 3) at distances of several kpc (Table 1), the linear spatial resolution is below 1000 AU for the most nearby sources (Table 3).

In contrast to likely not resolving all substructures within the regions, we nevertheless observe strong fragmentation in many targets. In particular, given the above estimated Jeans length between ~5500 and 27 700 AU (depending on density and temperature), most regions appear to fragment at or below this thermal Jeans length scale. Alternatively, the cores could have initially fragmented on Jeans length scales, and then the fragments could have approached each other even further due to the ongoing bulk motions from the global collapse of the regions. In contrast to that, the turbulent Jeans analysis, which includes the turbulent contributions to the sound speed, results in significantly larger mass and length scales (e.g., Pillai et al. 2011; Wang et al. 2014) than the classical thermal Jeans analysis.

Fig. 12 Nearest-neighbor separation histogram from minimum spanning tree analysis

6 Discussion

Fragmentation occurs in general on various spatial scales and is likely to be a hierarchical process. Within our CORE project, we investigated the fragmentation processes on clump scales in high-mass star-forming regions. We concentrated on the dense central structures on scales above ~1000 AU and roughly below 50 000 AU or 0.25 pc. These largest scales correspond roughly to the largest theoretically recoverable scales with 15 m baselines at 3 kpc distance (Sect. 4). In the continuum study presented here we investigate the fragmentation of clumps into cores. Fragmentation on smaller disk-like scales will also be investigated by the CORE program, however, this is more strongly based on the spectral line data and will be discussed in complementary papers (e.g., Ahmadi et al. 2018; Bosco et al. in prep.).

6.1 Thermal vs. turbulent fragmentation

With respect to the fragmentation of massive gas clumps, we identify some important questions: what controls the fragmentation properties of high-mass star-forming clumps? Is thermal Jeans fragmentation sufficient? How important are additional parameters like an initial non-uniform density profile or the magnetic field properties? How important is global accretion onto the clump from the diffuse interstellar medium (ISM)?

Regarding turbulent and thermal contributions, a number of studies have investigated this problem. For example, Wang et al. (2014) found that the observed masses of fragments within massive infrared dark cloud clumps are often more than 10 M_⊙. These masses are an order of magnitude larger than the thermal Jeans mass of the clump. Therefore, they argue that the massive cores in a protocluster are more consistent with turbulent Jeans fragmentation (i.e., including a turbulent contribution to the velocity dispersion). Similar results were found by Pillai et al. (2011) in their study of two young preprotocluster regions. On the other hand, Palau et al. (2013, 2014, 2015) found in their compiled sample of more evolved (IR-bright) star-forming regions that the masses of most of the fragments are comparable to the expected thermal Jeans mass, while the most massive fragments have masses a factor of 10 larger than the Jeans mass. Palau et al. (2013, 2014, 2015) concluded that these objects are consistent with thermal Jeans fragmentation of the parental cloud, in agreement with recent other investigations (e.g., Henshaw et al. 2017; Cyganowski et al. 2017). Recent ALMA studies of regions containing hypercompact HII regions also show small fragment separation scales (Klaassen et al. 2018). In addition to this, Fontani et al. (2016) argue that the magnetic field is important for the fragmentation of IRAS 16061–5048c1 (see also Commerçon et al. 2011; Peters et al. 2011).

In our sample of high-mass star-forming regions, including regions in an evolutionary stage comparable to those studied by Palau et al. (2013, 2014, 2015), we find that most of the fragment masses approximately agree with a plausible range of Jeans masses, and most nearest-neighbor separations are below the predicted scales of thermal Jeans fragmentation. To explore this in more detail, in Fig. 13 we plot the derived core masses against the nearest neighbor separation derived from the minimum spanning tree analysis. The solid and dashed lines show the relation between both thermal Jeans mass and Jeans length depending on density and temperature. In general, we do not find a clear trend between the two properties, and distance does not seem to be the primary factor in the observed scatter either. The figure also shows that for the plausible range of densities and temperatures (10⁵ to 10⁷ cm⁻³ and 10 to 100 K) the observed parameters are difficult to explain. One has to keep in mind that both observables are lower limits: the mass because of missing flux and the separation because of projection effects. Accounting for these effects, the measurements could shift a bit closer to the predicted lines, but could also shift sources parallel to them. For comparison, in the turbulent Jeans fragmentation picture, the sound speed is replaced by the velocity dispersion (e.g., Wang et al. 2014), which is typically a factor of five to ten higher than the thermal sound speed (see H₂ CO line width Δ v(H₂CO) in Table 3). Even if not all the observed line width is caused by pure turbulent motions, but also has contributions from organized motions due to, for example, large-scale infall, the regions clearly exhibit turbulent motions. Since the Jeans length and mass depend to the first and third power on the sound speed, respectively, replacing the thermal sound speed with the turbulent sound speed would shift the drawn correlations in Fig. 13 largely outside the observed box beyond the top-right corner. While we cannot conclude that thermal fragmentation explains everything, our data seem to refute that a turbulent contribution is needed if one applies a simple Jeans analysis for these spatial scales.

Several factors contributed to the apparent difference in fragmentation analysis between Wang et al. (2014) or Pillai et al. (2011) on the one side, and Palau et al. (2013, 2014, 2015) and the study here on the other side. First, the Wang et al. (2014) sample, incorporating data from Zhang et al. (2009) and Zhang & Wang (2011), has a typical 1σ mass sensitivity of 1 M_⊙. Therefore, lower mass fragments close to the global Jeans mass were not detected in these observations. Indeed, more sensitive observations from ALMA toward one of the objects in the sample, IRDC G28.34, revealed lower mass fragments (Zhang et al. 2015). Secondly, time evolution must play a role since fragmentation is a continuous process. As mentioned in Sect. 5.4, the separation scales between fragments may also change with evolutionary time. In the picture of globally collapsing clouds and gas clumps, one would expect larger fragment separation at early evolutionary stages. Then, during the ongoing collapse, the fragments may move closer together, following the overall gravitational contraction of the region. Therefore, the observed state of fragmentation only represents a snapshot in the time evolution. The less-evolved regions such as those in Wang et al. (2014) or Pillai et al. (2011) may present a deficit of low-mass fragments because the typical density of the cloud or clump is still lower so that a distributed low-mass protostar population may not have formed yet (e.g., Zhang et al. 2015). Furthermore, the more evolved objects such as those in this paper here have higher densities (Fig. 10), and therefore experience more fragmentation and are potentially more advanced in forming low-mass protostars.

In addition to the presented fragmentation properties, we point out that the nearest separations of cores are peaking around the spatial resolution limit of the observation (Fig. 12). Hence, fragmentation is also expected on even smaller scales. This can be investigated for this sample by higher spatial resolution observations with the future upgraded NOEMA (the baselines lengths are expected to be doubled), and for more southern sources with ALMA.

Recently, Csengeri et al. (2017) reported limited fragmentation for earlier evolutionary stages based on Atacama Compact Array data at 3.5″−4.6″. At the given spatial resolution and a mass sensitivity > 11M_⊙ they find that in 77% of their sample only three or fewer massive cores are found. However, because of the lower angular resolution and worse mass sensitivity, a direct comparison between their and this study is not possible. The data of Csengeri et al. (2017) are complemented with ALMA 12 m array data, and the combined dataset will be very valuable for comparison with the CORE project.

A different aspect to be considered is that the fragmentation properties likely change with spatial scale. Kainulainen et al. (2013, 2017) have shown for two filaments (the infrared dark cloud G11.11 and the Orion integral shape filament) that the fragmentation properties appear to show distinct signatures at different spatial scales. In particular for the infrared dark cloud, Kainulainen et al. (2013) argue that filament fragmentation dominates on large spatial scales (≥ 1 pc), whereas on smaller spatial scales thermal Jeans fragmentation takes over (~ 0.2 pc). With respect to our CORE sample, analyzing the filamentary properties on larger spatial scales is beyond the scope of this paper. However, it is clear that the CORE study deals with massive star-forming regions at high densities and not with the larger-scale, potentially filamentary clouds. In relation to the work by Kainulainen et al. (2013, 2017), our study is in the second regime that would be dominated by Jeans fragmentation. Therefore, our general result that the CORE sample is more consistent with thermal Jeans fragmentation is in agreement with the results by Kainulainen et al. (2013, 2017).

Fig. 13 Fragment masses against nearest neighbor separation from the minimum spanning tree analysis. The solid line corresponds to the Jeans lengths and Jeans masses calculated at 50 K for a density grid between 105 and 107 cm−3. For comparison, the dashed line corresponds to the Jeans lengths and Jeans masses calculated at a fixed density of 5 × 105 cm−3 (Beuther et al. 2002) with temperatures between 10 and 100 K. The color coding shows the distances of the sources.

6.2 Fragmentation diversity

Our sample clearly shows that the fragmentation properties within high-mass star-forming regions are not uniform, finding a diversity from highly fragmented regions to those that host one or only very few cores (see also, Bontemps et al. 2010; Palau et al. 2013; Csengeri et al. 2017). While this sample seems in general largely consistent with thermal Jeans fragmentation (see previous subsection), it should be noted that we also find a few massive cores in excess of 10 M_⊙ (Sect. 5.3 and Fig. 8). A high level of fragmentation with many low-mass cores favors high-mass star formation scenarios in the framework of competitive accretion (e.g., Clark & Bonnell 2006; Bonnell et al. 2007; Smith et al. 2009), whereas individual massive cores are more strongly needed in the turbulent core picture (e.g., McKee & Tan 2003; Tan et al. 2014). Because we find examples for both pictures in our CORE sample, this may indicate that different high-mass star formation scenarios are possible or even interplay with each other.

Since the sample is selected to host high-mass protostellar objects (HMPOs), the range of evolutionary stages is not broad. Nevertheless, as discussed in Sect. 2, within the HMPO category, we cover regions with varying IR-brightness and luminosity-to-mass ratios. Hence, while evolution is unlikely to be the main explanation for the observed fragmentation diversity, it cannot be entirely excluded. Furthermore, as discussed in the previous subsection, different levels of initial turbulence are also unlikely to be the underlying cause. Other possibilities to explain the different levels of fragmentation are variations in the initial density profiles and/or variations in the magnetic field properties. Differences in the density profiles could also arise from environmental effects like global collapse where the central gas clumps are continuously fed by some larger-scale cloud envelope.

Since the whole sample is observed with rather uniform uv-coverages, one wonders whether the amount of missing flux may be related to density structure of the parental gas clump and by that to the observed fragmentation properties of the cores. Therefore, we compare a few extreme cases: The two comparably isolated regions AFGL2591 and NGC 7538IRS1 (both at similar distances at 3.3 and 2.7 kpc, Table 1) show very different amounts of missing flux with values of 84% and 50% of the flux being filtered out. At the other extreme, two highly fragmented sources like S87IRS1 and IRAS21078 (at distances of 2.2 and 1.5 kpc, Table 1) also exhibit very different values of 87% and 53% of flux being filtered out. Hence, the overall fraction of flux being lost because of the interferometric observations – or rather the amount of mass in a diffuse, larger-scale reservoir – appears not to be an important issue for the observed fragmentation differences.

Girichidis et al. (2011) have used simulations of star-forming regions to show how the density profile affects the level of fragmentation: While flat profiles (ρ ∝ constant) resulted in many fragments, they find that density profiles like ρ ∝ r⁻² (over cloud radii of ~0.1 pc) quickly lead to the formation of a single object at the center where further fragmentation is prohibited. In their simulations, the intermediate case with ρ ∝ r^−1.5 is also dominated by a central object but additional fragments can form depending on their initial turbulence field. Observations of the density profiles of high-mass gas clumps by different groups typically find density slopes ρ ∝ r^−α with α between 1.5 and 2.6 (e.g., Beuther et al. 2002; Mueller et al. 2002; Fontani et al. 2002; Hatchell & van der Tak 2003). Furthermore, Palau et al. (2014) found a weak inverse trend between level of fragmentation and steepness of density profile, in other words, less fragmentation for steep density profiles. Hence, it seems reasonable that a range of initial density profiles can at least partly explain the observed diversity of fragmentation properties in our CORE sample. In future work, we will follow up on this and willinvestigate the density structure of the regions based on the combination of single-dish data with the interferometer data in more depth.

In addition to this, different magnetic field properties in the parental gas clumps can cause similar effects. Typically, the ratio between gravity and magnetic field is phrased in terms of the critical mass-to-flux ratio (e.g., Tilley & Pudritz 2007). Commerçon et al. (2011) modeled the collapse of high-mass star-forming regions with a range of magnetic field strengths. While their low-magnetic field case results in a larger number of fragments, the high-magnetic field case is dominated by a central massive object (see also Fontani et al. 2016). Similarly, Peters et al. (2011) also find reduced fragmentation with increasing magnetic field strength. To really differentiate whether the initial density profile and/or magnetic field properties are the dominant reason explaining the observed fragmentation diversity, we need to know the magnetic field strength, as well as the initial density profile. For two regions within the CORE sample (W3(H₂O) and NGC 7538IRS1) magnetic field studies have already been conducted with the Submillimeter Array on arcsecond resolution scale (Chen et al. 2012; Frau et al. 2014). The derived magnetic field strengths are comparably high in both regions with 17.0 and 2.5 mG, respectively. Since both regions exhibit very few or even only one fragment, the observed high magnetic field values are consistent with the low degree of fragmentation in these two regions. Future investigations in this direction are anticipated for the whole sample, which will reveal, in particular, whether regions with a high degree of fragmentation have a lower magnetic field strength.

7 Conclusions and summary

With the goals of studying the fragmentation, disk formation, outflows and chemical properties during the birth of the most massive stars, we have conducted the IRAM NOEMA large program CORE, observing a sample of 20 high-mass star-forming regions at 0.3″−0.4″ resolution in the 1.37 mm continuum and spectral line emission. In this paper, we present the survey scope, its main observational characteristics, the sample selection and the overall goals of the project. More details about the project as well as the first data release of the continuum data are provided online⁵. For a first scientific analysis of the data, we concentrated on the 1.37 mm dust continuum emission to investigate the fragmentation properties during early high-mass cluster formation.

We observed diverse fragmentation morphologies ranging from regions that are dominated by single high-mass cores to those that fragment into up to 20 cores. Since the sample contains mainly high-mass protostellar objects (although with some range of evolution within that category), larger-scale evolutionary effects are unlikely to explain all the differences. Observational artifacts like interferometric missing flux or different physical resolution can also be ruled out. The typical nearest neighbor separations peak below the thermal Jeans length determined from estimates of the initial average cloud density, indicating that thermal gravitational fragmentation is sufficient to explain the main observed core separations, and that additional turbulent contributions to the Jeans analysis are not needed for this sample. The diversity between regions with few or only one fragment vs. those with many fragments may be explained by differences in the initial density structures of the maternal gas clumps (potentially caused by environmental effects like global gas infall from a surrounding envelope) and/or variations in the initial magnetic field configurations. Since the nearest neighbor separation peaks around our spatial resolution limit, it is likely that further fragmentation takes place on even smaller spatial scales. With NOEMA, we will be able to address such questions for this northern hemisphere sample in a few years when the available baseline lengths will be doubled. Furthermore, ALMA observations of complementary southern hemisphere sources will investigate these questions in even greater depth.

Other scientific questions related to the disk formation, outflow properties and chemical processes during the formation of high-mass stars will be addressed by complementary CORE papers focusing on the spectral line data.

Acknowledgements

This work is based on observations carried out under project number L14AB with the IRAM NOEMA Interferometer and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This paper made use of information from the Red MSX Source survey database at http://rms.leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi which was constructed with support from the Science and Technology Facilities Council of the UK. H.B., A.A., J.C.M. and F.B. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648505. R.K. acknowledges financial support via the Emmy Noether Research Program funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1. R.G.M. acknowledges support from UNAM-PAPIIT program IA102817. T.C.s acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) via the SPP (priority program) 1573 “Physics of the ISM”. A.S.M. acknowledges support from Deutsche Forschungsgemeinschaft through grant SFB956 (subproject A6). A.P. acknowledges financial support from UNAM and CONACyT, México.

Appendix A Additional table

Appendix B Simulating CORE observations

To get a better quantitative understanding on how much the imaging and spatial filtering of the interferometer affects our results, we simulated CORE observations using a distance-scaled bolometer dataset from the Orion molecular cloud. The original Orion 850 μm data are thelarge-scale SCUBA-2 observations of Orion A by Lane et al. (2016).

For the simulations, several steps had to be applied to this original dataset. (i) Since our CORE data are at 1.3 mm wavelength, we scaled the flux densities S of the SCUBA-2 850 μm data by a typical spectral index at (sub)mm wavelength of S ∝ ν^3.5. (ii) We adapted the intrinsic spatial resolution of the data from the 450 pc distance of Orion (Lane et al. 2016) to 3 kpc (the typical distance of the CORE sample). This implied that the flux densities decrease by $S\propto {\left(\frac{450}{3000}\right)}^2$. (iii) The angular resolution of the original SCUBA-2 data of 14.6″ corresponds at the given distance of 450 pc to a linear resolution of 6570 AU. To maintain the same linear resolution at 3 kpc distance, we changed the pixel size of the image accordingly (from 3″ to 0.45″ pixel size). (iv) Furthermore, we converted the units of the image from mJy arcsec⁻² to mJy beam⁻¹ and then to K.

From the originally large Orion A image, we selected two subregions as examples: (a) a very strong region (Orion-KL), and (b) a fragmented weaker core within the northern part of the integral shape filament (refered to as “Orion-north” in the remaining part of the section). These two files are our model images that are run through the NOEMA simulator as described below. The two model images in units of mJy beam⁻¹ are shown in Figs. B.1 and B.2 (left panels).

Fig. B.1 Simulated CORE observations at 1.3 mm wavelength using original data obtained for Orion-KL by Lane et al. (2016). For details of the simulations see Appendix B. The beam for the left panel (Orion shifted to 3 kpc) is 2.19″, whereas the synthesized beam of the simulated image in the right panel is 0.44″ × 0.34″ (PA 44 deg). The contours in the left panel begin at 10 mJy beam−1 and continuein 30 mJy beam−1 steps. The contours in the right image start at 3σ and continuein 6σ steps (1σ ~ 0.3 mJy beam−1). The large circles are of 22″ diameter marking the FWHM of the primary beam.

Fig. B.2 Simulated CORE observations at 1.3 mm wavelength using original data obtained for a field in the northern filament of Orion A by Lane et al. (2016). For details of the simulations see Appendix B. The beam for the left panel (Orion shifted to 3 kpc) is 2.19″, whereas the synthesized beam of the simulated image in the right panel is 0.44″ × 0.34″ (PA 44 deg). The contours in the left panel are in 3 mJy beam−1 steps. The contours in the right image start are in 3σ steps (1σ ~ 0.04 mJy beam−1). The large circles are of 22″ diameter marking the FWHM of the primary beam.

The actual NOEMA simulations were then conducted within the GILDAS package. As a first step, a representative uv-coverage must be created within the ASTRO-subpackage by the command uv_tracks. To emulate our simulations best, we used a setup employing the D, B, and A-configurations. In each configuration the source was visited 15 times for a length of 15 min each with separations of 40 min between the visits. As target coordinates we used W3(H₂O) (Table 1). The resulting uv-coverage is shown in Fig. B.3.

Finally, the visibilities are produced within the GILDAS subpackage MAPPING and the task uv_fmodel using the above described model images (converted to K) and the given uv-coverage. The resulting visibility data files can then be imaged in exactly the same way as our original CORE datadescribed in Sect. 4. The synthesized beam of the simulated data is 0.44″× 0.34″ (PA = 44 deg).

As for the original CORE data, the rms varies corresponding to the side-lobes produced by the strongest sources in the field. In the strong Orion-KL field, the rms is 0.3 mJy beam⁻¹ whereas it is almost an order of magnitude lower (0.04 mJy beam⁻¹) in the weaker Orion-north field. For details of the simulation parameters see Table B.1.

Quantitatively, we extracted the fluxes toward the two main sources (#1 and #2) in Orion-north (Fig. B.2, Table B.1). Assuming optically thin dust emission at mm wavelengths with an assumed dust temperatureof 50 K, we were able to calculate the corresponding core masses as in the main part of the paper (Sect. 5). The approximate core masses of sources #1 and #2 in Orion-north are then 10.6 and 3.0 M_⊙, respectively. In the Orion-north field itself with the lower rms, we detect both sources well, however, with reduced integrated fluxes at 35%–40% ofthe model fluxes (Table B.1). Hence, the core masses are underestimated by the same factor. The situation would be different if these two sources were in the same field as Orion-KL. With an rms of 0.3 mJy beam⁻¹ in the Orion-KL simulation, the peak flux densities of sources #1 and #2 of 0.73 and 0.3 mJy beam⁻¹ are not above the 3σ thresholds.

Fig. B.3 Simulated uv-coverage.

If these Orion-north sources are typical for star-forming regions, then in fields with strong main sources, subsources of up ~10 M_⊙ could be difficult to detect. In contrast to that, in fields with weaker main sources, our core mass sensitivity should extend down to 1 M_⊙ or even lower. Therefore, our core-mass sensitivity within the CORE survey is dynamic-range and source-structure limited, and depends on the strongest source in the field. While for weaker sources like Orion-north, we reach almost the thermal noise limit, for stronger sources like Orion-KL, the dynamic range limit of these simulations is ~53.

It should be noted that the above derived core mass sensitivity strongly depends on the size of the cores and hence how much we resolve them. As shown in Table B.1, the point source mass sensitivity in both simulations is much smaller at 0.1 and 0.01 M_⊙, respectively. However, this is only valid for point sources. If our model cores were much more compact and all flux was within a single synthesized beam, we could easily detect them. But since we know from the original Orion observations their actual sizes, the emission extends over many beams which reduces the mass sensitivity for extended objects significantly.

If we now compare these simulations to the actual range of identified core masses in the survey (Fig. 8), we find that most cores have masses between ~0.2 and 6 M_⊙, below the regime estimated from our simulations. This difference can on the one hand be attributed to the missing flux implying that the core masses themselves are underestimated. But on the other hand, the source structure may also be different in our target regions compared to the Orion data used in these simulations. Beuther et al. (2015) have compared three different starless regions from the low-mass B68 to the intermediate-mass IRDC 19175 and the high-mass IRDC 18310-4 region. They found that the average densities between these three regions varied between 10⁴ cm⁻³ for the low-mass case to 10⁶ cm⁻³ for the high-mass region. While the Orion-north region within the integral-shape filament forms rather low- to intermediate-mass stars (in contrast to the high-mass region Orion-KL), our sample is selected to form high-mass stars. Hence, higher densities and more compact structures are expected forthe CORE sample. As outlined above, such compact structures make detections of lower-mass objects easier than if the emission is distributed over larger areas. Such structural differences may partly explain the differences between the simulated mass sensitivities and the observed ones.

Appendix C Individual continuum images

Fig. C.1 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.2 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.3 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.4 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.5 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.6 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

Fig. C.7 1.37mm continuum data for two CORE pilot sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.

References

All Tables

All Figures

	Fig. 1 Sample selection plot where the luminosity (in units of L⊙) is plotted against the MSX 21/8 μm color. Horizontal bars mark uncertainties in the color. While the blue sources fulfill our selection criteria, the red ones are below our luminosity cut of 104 L⊙. Green sources are those for which high-resolution mm data already exist and which were therefore excluded from the observations.
In the text

Fig. 2 Large-scale overview images for the whole CORE sample. The color scale shows three-color images with blue, green, and red from Spitzer 3.6, 4.5, and 8.0 μm for all sources except IRAS 23033, IRAS21078, G100, G094, and IRAS 23385 for which WISE 3.4, 4.6, and 12 μm data are presented. Furthermore, W3IRS4 uses Spitzer 3.6, 4.5 μm, and MSX 8 μm. The contours show SCUBA 850 μm continuum data (di Francesco et al. 2007; contour levels 20%, 40%, 60%, and 80% of the peak emission) for all sources except G100, G084, and G108 where these data do not exist.

In the text

	Fig. 3 Example wide-band spectrum extracted toward AFGL2591. The most important lines in the bandpass are marked.
In the text

	Fig. 4 Example uv-coverage for CepA. The different colors correspond to different observed (half-)tracks. Red and black correspond to D-array observations, blue and cyan to B-array, and green to A-array data.
In the text

	Fig. 5 Compilation of 1.37 mm continuum images for CORE sample on the same angular scale. The contouring is in 5σ steps (see Table 3). The sources are labeled in each panel, and the synthesized beams are shown at the bottom-left of each panel. A comparison figure converted to linear scales is shown in Fig. 6. Magnifications and absolute flux-scales are shown in Appendix C.
In the text

	Fig. 6 Compilation of 1.37 mm continuum images for CORE sample converted to linear resolution elements. The contouring isin 5σ steps (see Table 3). The sources are labeled in each panel, and the synthesized beams are shown at the bottom-left of each panel.
In the text

	Fig. 7 Number of identified cores plotted against the 5σ mass sensitivity. NGC 7538IRS1 is excluded because of its unusually high mass sensitivity limit.
In the text

	Fig. 8 Histograms of masses (top panel) and column densities (bottom panel) for all detected cores.
In the text

	Fig. 9 Gas column densities vs. masses for all detected cores. The color-coding shows the distances of the sources.
In the text

	Fig. 10 Mean core densities vs. masses for all detected cores. The color-coding shows the distances of the sources.
In the text

	Fig. 11 Core masses vs. equivalent radii for all detected cores. The color-coding shows the distances of the sources. The dashed lines show constant column density with levels of 1023, 1024, and 1025 cm−2 from right to left.
In the text

	Fig. 12 Nearest-neighbor separation histogram from minimum spanning tree analysis
In the text

Fig. 13 Fragment masses against nearest neighbor separation from the minimum spanning tree analysis. The solid line corresponds to the Jeans lengths and Jeans masses calculated at 50 K for a density grid between 105 and 107 cm−3. For comparison, the dashed line corresponds to the Jeans lengths and Jeans masses calculated at a fixed density of 5 × 105 cm−3 (Beuther et al. 2002) with temperatures between 10 and 100 K. The color coding shows the distances of the sources.

In the text

Fig. B.1 Simulated CORE observations at 1.3 mm wavelength using original data obtained for Orion-KL by Lane et al. (2016). For details of the simulations see Appendix B. The beam for the left panel (Orion shifted to 3 kpc) is 2.19″, whereas the synthesized beam of the simulated image in the right panel is 0.44″ × 0.34″ (PA 44 deg). The contours in the left panel begin at 10 mJy beam−1 and continuein 30 mJy beam−1 steps. The contours in the right image start at 3σ and continuein 6σ steps (1σ ~ 0.3 mJy beam−1). The large circles are of 22″ diameter marking the FWHM of the primary beam.

In the text

Fig. B.2 Simulated CORE observations at 1.3 mm wavelength using original data obtained for a field in the northern filament of Orion A by Lane et al. (2016). For details of the simulations see Appendix B. The beam for the left panel (Orion shifted to 3 kpc) is 2.19″, whereas the synthesized beam of the simulated image in the right panel is 0.44″ × 0.34″ (PA 44 deg). The contours in the left panel are in 3 mJy beam−1 steps. The contours in the right image start are in 3σ steps (1σ ~ 0.04 mJy beam−1). The large circles are of 22″ diameter marking the FWHM of the primary beam.

In the text

	Fig. B.3 Simulated uv-coverage.
In the text

	Fig. C.1 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.2 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.3 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.4 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.5 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.6 1.37mm continuum data for CORE sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text

	Fig. C.7 1.37mm continuum data for two CORE pilot sources. The left panels always show the data without applying self-calibration, and the right panels show them after applying self-calibration. The contours are always in 5σ steps (see Table 3). The right panels indicate the cores identified with clumpfind.
In the text