Fragmentation and kinematics in high-mass star formation: CORE-extension targeting two very young high-mass star-forming regions

Context. The formation of high-mass star-forming regions from their parental gas cloud and the subsequent fragmentation processes lie at the heart of star formation research. Aims. We aim to study the dynamical and fragmentation properties at very early evolutionary stages of high-mass star formation. Methods. Employing the NOrthern Extended Millimeter Array (NOEMA) and the IRAM 30m telescope, we observed two young high-mass star-forming regions, ISOSS22478 and ISOSS23053, in the 1.3mm continuum and spectral line emission at a high angular resolution ( ∼ 0.8 (cid:48)(cid:48) ). Results. We resolved 29 cores that are mostly located along ﬁlament-like structures. Depending on the temperature assumption, these cores follow a mass-size relation of approximately M ∝ r 2 . 0 ± 0 . 3 , corresponding to constant mean column densities. However, with di ﬀ erent temperature assumptions, a steeper mass-size relation up to M ∝ r 3 . 0 ± 0 . 2 , which would be more likely to correspond to constant mean volume densities, cannot be ruled out. The correlation of the core masses with their nearest neighbor separations is consistent with thermal Jeans fragmentation. We found hardly any core separations at the spatial resolution limit, indicating that the data resolve the large-scale fragmentation well. Although the kinematics of the two regions appear very di ﬀ erent at ﬁrst sight – multiple velocity components along ﬁlaments in ISOSS22478 versus a steep velocity gradient of more than 50kms − 1 pc − 1 in ISOSS23053 – the ﬁndings can all be explained within the framework of a dynamical cloud collapse scenario. Conclusions. While our data are consistent with a dynamical cloud collapse scenario and subsequent thermal Jeans fragmentation, the importance of additional environmental properties, such as the magnetization of the gas or external shocks triggering converging gas ﬂows, is nonetheless not as well constrained and would require future investigation.


Introduction
During the hierarchical process of star formation, the gas has to flow from the large cloud scales down to the smallest scales, where disks and protostars form. While questions regarding the small-scale structure and the presence of disks toward more evolved high-mass protostellar objects have been subject to intense study for several decades (e.g., Cesaroni et al. 2007;Beltrán & de Wit 2016;Beuther et al. 2018), the earliest evolutionary stages of collapse and fragmentation during the formation processes of high-mass stars have not been investigated in depth thus far. This is in part related to the short evolutionary time-scales of early massive star formation (e.g., Beuther et al. 2007;Zinnecker & Yorke 2007;Tan et al. 2014) and the fact that hardly any high-mass starless core candidates manage to survive further scrutiny when additional star formation indicators become detected (e.g., Motte et al. 2007;Tan et al. 2016;Feng et al. 2016). A few exceptional regions exist that have remained starless even with high-resolution millimeter interferometric observations (e.g., G11.92 Cyganowski et al. 2014Cyganowski et al. , 2017IRDC18310-4 Beuther et al. 2013, 2015aW43-MM1 Nony et al. 2018). However, at least the latter two examples are at comparably large distances of 3.4 and 4.9 kpc, respectively, reducing the ability to investigate structures at small physical spatial scales. Furthermore, recent investigations with the Atacama Large Millimeter Array (ALMA) and the Submillimeter Array (SMA) also started to further dissect samples of very young high-mass star-forming regions (e.g., Sanhueza et al. 2019;Svoboda et al. 2019;Li et al. 2019Li et al. , 2020a. Important questions regarding the earliest evolutionary stages of high-mass star formation are related to the fragmentation properties of the dense gas and the associated kinematic properties. For example, we consider whether these very young regions fragment early on into many lower-mass cores that may then competitively accrete gas from the surrounding envelope (e.g., Bonnell et al. 2007;Smith et al. 2009). Alternatively, they may fragment into more massive cores with potentially less lowmass cores early on (e.g., Tan et al. 2014;Zhang et al. 2015;Csengeri et al. 2017;Beuther et al. 2018)? We ask how important filamentary accretion processes are (e.g., Schneider et al. 2010;Peretto et al. 2014;André et al. 2014;Chira et al. 2018;Hennebelle 2018;Padoan et al. 2020) and what the relative importance is of cloud-scale gravo-turbulent or thermal fragmentation versus fragmentation of the gravitationally unstable accretion flows for the star formation process (e.g., Peters et al. 2010)? Furthermore, it has been argued that high-mass star formation starts and proceeds in turbulent gas clumps (e.g., McKee & Tan 2003). However, recent high-spatial-resolution studies of a few very young high-mass star-forming have regions revealed that at high enough spatial and spectral resolution, the spectral lines re-A&A proofs: manuscript no. 40106  (Ragan et al. 2012;Di Francesco et al. 2008). The contour levels are from 18 to 98% of the peak emission of 1.88 Jy beam −1 and 14.88 Jy beam −1 for the two regions, respectively. Linear scale bars are shown and the red circles outline the NOEMA mosaics with a primary beam FWHM of 22 . solve into multiple components, each with line widths that are nearly consistent with thermal line broadening (e.g., Beuther et al. 2015a;Hacar et al. 2018;Li et al. 2020b). Therefore, it may well be the case that the individual cores do not have turbulent properties but it is, rather, only the superposition of many of these individual cores and spectral components mimicking the broader lines that are typically seen with single-dish telescopes at lower spatial resolution (e.g., Smith et al. 2013;Hacar et al. 2017b). At a somewhat coarser resolution of >10 000 au, similar results indicating sub-sonic non-thermal motions were also found by Ragan et al. (2015) and Sokolov et al. (2018) in lowerdensity lines of N 2 H + and NH 3 , respectively.
With the goal of studying the kinematic and fragmentation properties at the onset of high-mass star formation, we employ high spatial resolution and high-density tracer observations at mm wavelengths toward two regions at very early evolutionary stages. The targets are part of the IRAM NOEMA large program CORE aimed at studying the fragmentation, disk properties, kinematics, and chemistry during the formation of highmass stars (Beuther et al. 2018). In the following, we summarize the main aspects of the CORE program.

The CORE and CORE-extension projects
With the goal to study the fragmentation, disk, outflow, and chemical properties of high-mass star formation, we have embarked on the NOEMA large program CORE, which observes 20 high-mass star-forming regions at a high spatial resolution (0.3 − 0.4 ) in the 1.3 mm line and continuum emission. The sample was selected based on being mainly in the evolutionary stage of high-mass protostellar objects (HMPOs), also known as MYSOs (massive young stellar objects). An overview and more details about the program can be found in Beuther et al. (2018) 1 .
To broaden the scope of CORE to even younger evolutionary stages, we conducted the CORE-extension part of the project, observing two much younger regions (ISOSS22478 and ISOSS23053) also at high spatial resolution and with a similar 1 See also http://www.mpia.de/core spectral setup. The goal of this CORE-extension is to investigate earlier evolutionary stages with a particular focus on the gas dynamics and chemical properties at the onset of high-mass star formation. Selection criteria for their youth were mainly their low temperatures and luminosities as well as the low luminosity-to-mass ratios L/M (∼0.4 and ∼2.2 for ISOSS22478 and ISOSS2305, respectively). These L/M ratios are typically more than an order of magnitude below those of the original CORE sample (Beuther et al. 2018). Below, we introduce the target regions and Sect. 2.2 presents the technical details of the CORE-extension project.

Target regions
While the main part of the program targets 20 HMPOs with already embedded massive young stellar objects, we also observed two younger regions in slightly more extended mosaics to investigate the earliest formation processes. These two regions were originally identified via the ISO Serendipity Survey (ISOSS), which observed the sky serendipitously also during slew times with the Infrared Space Observatory (ISO) at 170 µm (Bogun et al. 1996;Krause et al. 2004). Since the regions are located in the outer Galaxy, and not in front of strong Milky Way background emission, they do not appear as infrared dark clouds (IRDCs). However, their general physical properties put them in a similarly young evolutionary stage as typical IRDCs.
ISOSS22478+6357: Earlier studies of this region were conducted by Hennemann et al. (2008), Ragan et al. (2012) and Bihr et al. (2015). At a distance of ∼3.23 kpc, the region has a mass reservoir of ∼140 M , placing it in the category of intermediatemass star-forming regions. The average temperatures based on NH 3 and dust observations are between ∼13 K and ∼21 K, and the region still has a low total luminosity of ∼55 L (Ragan et al. 2012;Bihr et al. 2015). Figure 1 (left panel) gives a large-scale overview of the region comprising the Herschel 70 µm emission and the SCUBA 450 µm data (Ragan et al. 2012;Di Francesco et al. 2008).
ISOSS23053+5953: This region has also previously been investigated, for example, by Birkmann et al. (2007), Ragan et al. (2012) and Bihr et al. (2015). At a distance of ∼4.31 kpc, the region has a mass reservoir and luminosity of 610 M and 1313 L , respectively. Average NH 3 and dust temperature for the region are ∼18 K and ∼22 K (Bihr et al. 2015;Ragan et al. 2012). Figure  1 (right panel) shows the corresponding large-scale overview image of the region. While early studies already showed signatures of infall (Birkmann et al. 2007), an analysis of the kinematics of this region based on VLA (Very Large Array) NH 3 observations revealed indications of two velocity components that may collide at the positions of the central cores (Bihr et al. 2015).

The CORE-extension technical setup
While the original CORE program utilized the old receiver and backend system with only a ∼4 GHz wide band in a single sideband between 217.167 and 220.834 GHz, the upgraded NOEMA now allows us to use, in dual-polarization, a much broader bandwidth in the double-sideband mode with 7.8 GHz in each sideband. Specifically, we are covering the frequency ranges from ∼213.3 to ∼221.1 GHz as well as from ∼228.8 to ∼236.6 GHz. This spectral setup covers a multitude of strong spectral lines from many important molecules, for instance, deuterated species such as DCO + , DNC, or N 2 D + , the temperature tracers H 2 CO or CH 3 CN, CO, and its isotopologues 13 CO and C 18 O, shock tracers like SiO or larger molecules such as CH 3 OH or HCCCH. While the entire bandpass is covered at low spectral resolution (2 MHz corresponding to ∼2.7 km s −1 at 220 GHz), all important lines are also covered by separate high-spectral-resolution units (0.06 MHz corresponding to ∼0.08 km s −1 at 220 GHz).
More details about the spectral setup and the chemical properties of the regions will be given in a future paper by Gieser et al. (in prep.). Here, we focus on the dynamical properties and potential gas flows in these two regions.

Observations and data
The two regions were observed with NOEMA in February and March 2019 with ten antennas in the A, C, and D configurations covering baseline lengths roughly between 18 and 774 m. The sources ISOSS22478 and ISOSS23053 were observed as six-field and four-field mosaics, respectively (see red circles in Fig. 1 Calibration and imaging of the data was performed with the clic and mapping software of the gildas package 2 . To achieve optimal imaging quality for these mosaics, we applied natural weighting during the imaging process. To create the continuum data, only the line-free parts of the spectrum were used. Since the regions are comparably line-poor, only small parts of the entire bandpass (Sect. 2.2) had to be excluded. Specifically, we excluded 280 MHz and 1800 MHz of the entire bandpass of ∼15 600 MHz for ISOSS22478 and ISOSS23053, respectively. For the 1.3 mm continuum data, this resulted in synthesized beams of 0.92 × 0.73 (P.A. 51 deg) and 0.84 × 0.74 (P.A. 67 deg) for ISOSS22478 and ISOSS23053, respectively. This corresponds to approximate linear resolution elements of ∼2700 and ∼3500 au, respectively. The corresponding 1σ continuum rms values are 0.057 mJy beam −1 and 0.16 mJy beam −1 , respectively. The absolute flux scale is estimated to be correct within 20%.
For the spectral line data, we obtained complementary single-dish observations with the IRAM 30 m telescope to compensate for the missing short spacings. These 30 m observations were conducted in the on-the-fly mode typically achieving rms values of ∼0.1 K (T mb ).
The NOEMA and 30 m data were then combined in the imaging process with the task uv_short, and we present the merged data for DCO + (3-2), SiO  and several H 2 CO lines at a spectral resolution of 0.3 km s −1 . Again, natural weighting was applied, and the final 1σ rms and synthesized beam values are presented in Table 1. The other spectral lines and a detailed chemical analysis of the regions will be presented in Gieser et al. (in prep.).

Fragmentation and continuum emission
To investigate the fragmentation properties of the large-scale gas clumps visible in Figure 1, we show the 1.3 mm dust continuum emission of the NOEMA data at a spatial resolution of ∼0.8 (Table 1) in Figure 2. Both regions clearly fragment into numerous cores. Because we have no short spacing information for the continuum emission, these cores appear to be separated from each other. As we go on to see in the spectral line data discussed below (Sect. 4.2), these cores are nevertheless connected by filamentary gas structures. Furthermore, we find that all cores are associated with 70 µm emission, hence, none of them are genuinely starless any longer. Although the number of cores is not excessively large, to systematically derive the flux densities, offsets and sizes of the cores, we use the clumpfind algorithm originally introduced by Williams et al. (1994). As a low-intensity threshold, we used the 5σ contours presented in Fig. 3. We identify 15 and 14 cores for ISOSS22478 and ISOSS23053, respectively. Peak flux densities S peak , integrated fluxes S , R.A. and Dec., offset positions (∆x, ∆y) and equivalent core radii are presented in Table 2. The equivalent core radii are calculated from the measured core area assuming a spherical distribution. No deconvolution from the synthesized beam was applied.
Core masses and peak column densities can be estimated from the integrated fluxes, S , and the peak flux densities, S peak , assuming optically thin dust continuum emission. Following Hildebrand (1983) and Schuller et al. (2009), we use a gas-todust mass ratio of 150 (Draine 2011) and a dust mass absorption coefficient κ = 0.9 cm 2 g −1 (Ossenkopf & Henning 1994 at densities of 10 6 cm −3 with thin ice mantles).
The last bits of information needed is the temperature of the gas and dust. For both regions several temperature estimates exist. Ragan et al. (2012) used Herschel far-infrared data and fitted the spectral energy distributions. This resulted in average dust temperatures for ISOSS22478 and ISOSS23053 of ∼21 and ∼22 K, respectively. Both regions were also imaged with the Very Large Array in NH 3 emission, and average temperatures from these data are 13 and 18 K for the two regions, respectively (Bihr et al. 2015). Within this new NOEMA dataset we also have the H 2 CO lines around 218 GHz that can be used as thermometer as well (Mangum & Wootten 1993). Details about the temperature derivation and structure are presented in Sect. 4.2. Relevant for our mass and column density analysis here is that for ISOSS22478 the temperatures range approximately from values around 10 K to roughly 50 K (Table 2). For ISOSS23053, the range is broader, from around 10 K to more than 100 K. In Sect. 4.2, we discuss in more detail that a lot of the high temperature gas also appears to be associated with shocks.
The above three temperature estimates indeed show some variance. Since our mass and column density estimates are based on dust continuum observations of dense gas cores, the Herschel dust, and the NH 3 estimates may give good constraints. How-H. Beuther et al.: Kinematics of cloud collapse during high-mass star formation 15 Notes: * The nominally derived temperature was either slightly below 10 K or could not be well fitted (#10). We set a 10 K lower floor here. ever, the new NOEMA H 2 CO temperature data have a higher spatial resolution and are also able to trace a higher dynamic range in temperature than the dust and NH 3 (1,1)/(2,2) data. As discussed in Sect. 4.2, the H 2 CO temperatures appear to be also affected by the shock structure. While the dust and NH 3 data may underestimate the temperatures toward the main cores, the H 2 CO data may overestimate temperatures toward lower-mass cores in the regions that are close to shock enhancements of the temperatures (e.g., core #12 in ISOSS22478 or core #14 in ISOSS23053, Table 2). Hence, it is not a priori clear which of the temperature estimates are better for determining the masses and column densities. Therefore, in the following, we derive our mass and column densities estimates with two assumptions. Based on the dust and NH 3 temperature measurements, we use uniformly 20 K for all mass and column density estimates in the first approach. In the second, we use the temperatures derived from the H 2 CO data at each core position for the mass and column density estimates. We refer to the two approaches as T 20K and T H2CO . While the interferometrically filtered out emission in the continuum data implies that the measured fluxes are lower limits, considering the additional uncertainties in the temperature estimate, distances, dust absorption coefficient κ as well as the gas-to-dust mass ra-tio, we estimate the masses and column densities to be accurate within a factor of 2-4 (see also the CORE study of the original sample, Beuther et al. 2018). The derived core masses (M 20K and M T H2CO ), peak column densities ( N peak,20K and N peak,T H2CO ), and H 2 CO temperatures (T H2CO ) for the two approaches outlined above are presented in Table 2 for all 29 cores (15 and 14 cores for ISOSS22478 and ISOSS23053, respectively). The range of core masses for ISOSS22478 is between ∼0.05 and ∼4.5 M whereas core masses for ISOSS23053 range between ∼0.3 and 39.5 M , depending also on the temperature assumptions. While the difference in ranges may partially be attributed to the slightly larger distance of ISOSS23053, this region is also far more massive in general (Sect. 2.1). Hence, it is not too surprising that more massive cores are detected in that region. The column densities for both regions range between a few times 10 22 to ∼10 24 cm −2 . With the given assumptions at 20 K the 3σ mass and column densities for ISOSS22478 are 0.06 M and 1.5 × 10 22 cm −2 , respectively. The corresponding 3σ values again at 20 K for ISOSS23053 are 0.3 M and 4.1 × 10 22 cm −2 .
The total mass in all cores for ISOSS22478 and ISOSS23053 are 14.7 and 121.9 M assuming a temperature of 20 K, respec- Fig. 4. Relations of core parameters. The top panels show the peak column density and mean density plotted against the mass while the bottom panels present the mean density and mass plotted against the equivalent radius. The red and blue colors always show the data for ISOSS22478 and ISOSS23053, respectively. The full pentagons show the results estimated assuming constant temperatures (T 20K ) whereas the crosses present the data assuming the temperatures derived from the H 2 CO line data (T H2CO ). The bottom-left panel presents as green line a fit to the data following the T H2CO approach. The bottom-right panel also presents fits to the data (green continuous lines) and lines of constant mean column density and constant mean density (dashed and dotted lines, respectively) as labeled in the plot. Respective error-bars in each panel correspond to a factor of 2 uncertainty for masses, column densities, and mean densities, as well as 1 kpc distance uncertainty for the equivalent radii (resulting in the different error-bar end markers for the radii).
tively. Comparing these values with the total mass reservoir available in the regions (Sect. 2.1), because there is no shortspacing information for the continuum data, at most 10% and 20% of the continuum flux can be recovered by the NOEMA observations for ISOSS22478 and ISOSS23053, respectively. Figure 4 presents several relations derived for the cores of the two regions. All plots show the derived parameters with both temperature assumptions. The top-left panel shows the column density against the core mass, and we find a clear correlation between the two parameters, independent of the temperature assumption. While for the T H2CO approach, the masses and column densities between the two regions overlap strongly, in the T 20K approach, the estimated masses and column densities are larger for ISOSS23053. Nevertheless, also in the T 20K approach, there is clear overlap between both regions and the general trend for these two quantities agrees well.
Employing the equivalent radius of the cores from the clumpfind analysis, we can also derive mean densities for all cores assuming spherical symmetry. The top-right and bottomleft panels of Fig. 4 plot these mean densities versus the core mass and the core size, respectively. While the scatter in the T 20K approach for the densities is smaller around almost constant H. Beuther et al.: Kinematics of cloud collapse during high-mass star formation densities of ∼ 10 6 cm −3 , the bottom-left panel of Fig. 4 is indicative of a decreasing trend of densities with increasing equivalent core radius in the T H2CO approach. Fitting a power-law to the size-density relation in this T H2CO approach, we find a relation of n ∝ r −1.0±0.3 . We note that in the classical 3rd Larson relation, Larson (1981) found that for his sample of CO molecular clouds the mean density is inversely related to the size (n ∝ r −1.1 ). The bottom-right plot of Fig. 4 shows the core masses against the equivalent radius, that is, the mass-size relation. While we see a good correlation between these two quantities, the slopes for the relation derived under the two temperature assumptions varies significantly. The mass-size relation in the T H2CO approach follows M ∝ r 2.0±0.3 , close to the lines of constant column densities as would be expected from Larson's relations (e.g., Heyer et al. 2009;Lombardi et al. 2010;Ballesteros-Paredes et al. 2019. In contrast to this, the mass-size relation in the T 20K approach follows a steeper relation of M ∝ r 3.0±0.2 , which is close to the line of constant density at 10 6 cm −3 . If the two regions are separated in the T 20K approach, the mass-size relations for ISOSS22478 and ISOSS23053 have power-law dependences of r 2.4±0.2 and r 2.6±0.2 , respectively. We return to this point in Sect. 5.1. To investigate the projected separations between cores, we employed a minimum spanning tree analysis using the astroML software package (VanderPlas et al. 2012) that determines the shortest possible distances connecting each of the cores. The derived nearest-neighbor separations are shown in Fig. 5 for both regions combined, and the basic results are listed in Table 3. While projection causes the observed nearest-neighbor separations to be at their lower limits, as compared to their real separation in three dimensions, the sensitivity, and spatial resolution limits that the observation are subject to, implying that the projected nearest-neighbor separations are, rather, likely to be the upper limits. In reality, having more cores should qualitatively reduce the real projected nearest-neighbor separations. Nevertheless, with a linear spatial resolution of approximately 2700 and 3500 au, it is interesting to note that barely any separations are found at our resolution limit. This is different in comparison to the original CORE sample of more evolved HMPOs where the nearest-neighbor histogram shows a strong peak near the resolution limit. While also for the younger ISOSS sources, we would expect that at higher spatial resolution a few of the cores are likely to split into multiple systems and, hence, some smaller separations would occur, the current data nevertheless show that we are resolving the large-scale fragmentation of the parental gas clumps well. A decrease in spatial separations during the evolution of a star-forming cluster is expected during the collapse of the gas clump. During the collapse, the densities increase and in the classical Jeans picture, the Jeans length and the corresponding separations decrease. In that scenario, the generally larger separations of the younger ISOSS sources studied here can naturally be explained.
It is interesting to note that in a study of the integral shape filament in Orion, Kainulainen et al. (2017) found that significant grouping of dense cores is found below scales of 17000 au (see also Román-Zúñiga et al. 2019), similar to what we find  2017) infer a second increase in source grouping below separations of 6000 au, which is not visible in our data. However, smaller-scale grouping (below our resolution limit here) was found in the original CORE data at approximately twice the spatial resolution with a peak at roughly 2000 AU (Beuther et al. 2018). This missing smallerscale grouping in the two regions studied here may potentially be attributed to insufficient spatial resolution. We come back to that in Sect. 5.1.

Kinematics and temperatures of the gas
We consider what the kinematics of the gas tell us about the formation processes within these two star-forming regions. To look at the dynamics in more detail, in the following we focus on the DCO + (3-2) emission that traces the gas properties at the given early evolutionary stage extremely well (e.g., Gerner et al. 2015). Figures 6 and 10 present an overview of the DCO + data. The left panels show the integrated intensities and the middle and right panels present the 1st and 2nd moment maps (intensity-weighted peak velocities and velocity dispersions), respectively. These are always the merged datasets combining the NOEMA ACD configurations with the IRAM 30 m data (Sect. 3).

ISOSS22478
For ISOSS22478, the integrated emission clearly shows extended filamentary gas emission structures that connect the 1.3 mm continuum cores. The 1st moment map in the middle panel of Fig. 6 reveals a complex velocity structure with slightly different velocity components for the different sub-regions. We return to this point below. The velocity dispersion shown in the right panel shows on average low values largely around or even below 1 km s −1 (see also the cases of G29.96 and G35.20 for similarly narrow lines in deuterated ammonia, NH 2 D, Pillai et al. 2011). For comparison, Fig. 7 presents individual example spectra toward six cores. While the spectrum toward core #2 is better fitted with two components, some other spectra are still consistent with single-component fits. The full width half maximum    7. DCO + (3-2) example spectra toward the labeled core positions in ISOSS22478 (Fig. 3). The spectra are shifted on the Y-axis for better presentation. The spectra of core #5 and core #6 are multiplied by 3 and 2 for clarity. The green lines show Gaussian fits, the FWHM linewidth ∆ are presented for all spectra. For core #2, two components were fitted.
A different way to look at the kinematics of the gas is to consider them as position-velocity slices along specific cuts within a region. Here, we concentrate on three slices marked in the middle panel of Fig. 6. Slice 1, connecting mainly cores #2 to #1 and #11, exhibits interesting multiple velocity components towards most of the core positions along the slice. Similarly, slice 3 also shows several velocity components at individual positions along the cut. In contrast to these two slices, the almost perpendicular cut, connecting mainly cores #6, #11, and #5 shows a different velocity gradient, with initially decreasing and then increasing Another important parameter is the temperature of the gas in the region. Figure 9 presents the temperature map we derived from the H 2 CO lines around 218 GHz (Table 1), which are a well-characterized thermometer of the interstellar gas (e.g., Mangum & Wootten 1993;Rodón et al. 2012;Gieser et al. 2019). To derive this temperature map, we used the eXtended CASA Line Analysis Software Suite (xclass, Möller et al. 2017 3 ) tool within the casa software package (McMullin et al. 2007). In xclass the spectra are fitted pixel by pixel taking also the optical depth into account. The mean uncertainties for the temperatures are on the order of 30%. More details on the xclass fitting are given in Gieser at al. (subm. to A& A). In addition to this, an in-depth analysis of the chemical and temperature structure of cores within the the two regions will be presented in Gieser et al. (in prep.). As we can see in Fig. 9, the temperatures in ISOSS22478 are homogeneously very low, typically below 20 K, and only in the central part of the map do they rise above 20 K over some area. The temperature enhancement above 50 K in the small condensation in the south of the map (near the green contours in Fig. 9) is considered to be real because it is the only spot where also SiO(5-4) emission is detected. Hence, while we find at one localized position a temperature enhancement that seems to be associated with shocks, as suggested by the detection of SiO emission, we do not identify strong temperature differences along the main filamentary parts of the cloud. This is set into contrast to the SiO emission and temperature structure in ISOSS23053 below.

ISOSS23053
In comparison to ISOSS22478, the DCO + gas emission in ISOSS23053 looks very different (Fig. 10). Visually, a filamentary structure is less apparent. In comparison to ISOSS22478, where the DCO + and mm continuum emission spatially agree relatively well, for ISOSS23053 many of the DCO + peaks are offset from the main mm continuum peak positions. At the given spatial resolution, the corresponding continuum peak brightness temperatures are all below 1 K and the dust continuum emission is optically thin. Hence, dust optical depth effects cannot explain these offsets. More important in the context of the dynamical state of the region, the 1st moment map in Fig. 10 exhibits two very distinct velocity components that merge in the middle of the region on scales below 0.1 pc, right where we find most of the strong mm continuum cores. These two velocity components were already identified in previous NH 3 data, and Bihr et al. (2015) suggested that they may be the signature of two converging gas streams that collide and trigger a new generation of star formation.
In that context, it is interesting to investigate the velocity dispersion. The 2nd moment map in Fig. 10 reveals that the velocity dispersion is largely low, at or below 1 km s −1 . However, we see an increase in velocity dispersion above 1 km s −1 towards the main cores, and even higher values a bit to the north of the two main central cores (#2 & #3). Comparing these moment maps to a few selected spectra (Fig. 11), most of the spectra have rather narrow full width half maximum ∆ values below 1.5 km s −1 with the exception of cores #1 and #2 that reveal ∆ above 2 km s −1 . The spectra also show clearly the two velocity components, for instance, core #7 peaks around -54 km s −1 whereas the emission peak of core #8 is rather at -51 km s −1 .
Regarding the two velocity components in ISOSS23053, Fig. 12 shows two position-velocity slices with the directions outlined in the middle panel of Fig. 10. Both slices clearly show the two velocity components that connect on very small spatial scales. While Bihr et al. (2015) report this velocity step as unresolved in their ∼4 NH 3 observations and give a lower limit to the gradient of > 30 km s −1 pc −1 , our higher spatial resolution allows us to better quantify this structure. For slices 1 and 2, we measure velocity gradients of ∼48 km s −1 pc −1 and ∼54 km s −1 pc −1 , respectively. They can be even larger in the SiO emission discussed in the following.
In considering shocks by potential converging gas flows, we also investigate the SiO(5-4) emission presented in Fig. 13. Here, SiO is expected to be strongly enhanced by shocks because Si gets sputtered from the grains and quickly forms SiO in the gas phase (e.g., Schilke et al. 1997;Anderl et al. 2013). Interestingly, the SiO emission is comparably weak toward the main cores in the region, but it shows enhanced emission at several positions throughout the cloud. The strongest integrated SiO emission is identified between the main cores #2&#3 and #6. Additional significant SiO emission is found between the cores #1 and #5 as well as in some elongated structures toward the north of the region (e.g., one labeled as "stream" in Fig. 13). While the peak velocities in the 1st moment map show the same bluered dichotomy as the previously shown by the DCO + data, the SiO linewidth increase visible in the 2nd moment map in Fig. 13 is very pronounced. We find increased velocity dispersion again A&A proofs: manuscript no. 40106 Fig. 10. NOEMA+30 m DCO + (3-2) data towards ISOSS23053. The three panels show in color scale the integrated intensity, the 1st and 2nd moment maps (intensity-weighted peak velocities and velocity dispersions), respectively. These maps were produced by clipping all data below an approximate 3σ threshold of 15 mJy beam −1 . The contours show the 1.3 mm continuum emission in 3σ steps (1σ ∼ 0.16 mJy beam −1 ) from 3 to 15σ. The beam is shown in the bottom-right of all panels, and a linear scale bar is shown in the middle panel. The arrows in the middle panel show the two position velocity slices presented in Fig. 12. To reduce noise signatures at the edges of the mosaic, for the 2nd moment map, we blanked the mosaic edges. The cores are labeled in the left panel. Fig. 11. DCO + (3-2) example spectra toward the labeled core positions in ISOSS23053 (Fig. 3). The spectra are shifted on the Y-axis for better presentation. The spectra of core #1, #2, and #7 and are multiplied by 2, 2, and 3 for clarity. The green lines show Gaussian fits and the FWHM linewidth ∆ is presented for all spectral lines. between cores #2/#3 and #6 as well as between the cores #1 and #5 and the northern stream.
While SiO is often associated with shocks from outflows (e.g., Chandler & Richer 2001;Palau et al. 2006;Cabrit et al. 2012), shocks in a more general sense also from other drivers are sufficient (e.g., Anderl et al. 2013). For example, shocks during cloud formation may contribute to the SiO formation and emission as well (e.g., Jiménez-Serra et al. 2010. In the case of ISOSS23053, the stream-like features exhibit elongated structures that are most likely caused by outflows. The orientation of the stream-feature labeled in Fig. 13 is towards core #3, which may be the driver of that outflow structure. However, other strong SiO emission peaks, in particular those between the cores #1 and #5 as well as between cores #2&#3 and #6, do not necessarily need to be of protostellar outflow origin. Although core #1 is reported to drive an outflow (Birkmann et al. 2007), the strong SiO emission is rather centered between cores #1 and #5. All the SiO features between the cores #1/#5 and #2/#3/#6 are spatially closely linked to the velocity jump best visible in Fig. 10 or nearby core #6 in the SiO emission in Fig. 13. Conducting a position velocity cut approximately through core #6 (middle panel of Fig. 13), the corresponding position-velocity diagram is presented in Fig. 14. The SiO emission along this cut covers a much broader velocity range than the DCO + emission shown in Fig. 12. Furthermore, the two velocity components appear even more clearly separated in space.
While it cannot be entirely excluded that these velocity differences are of protostellar outflow origin, the clear spatial separation of the blue-and red-shifted emission without a clear driving source at its center make the outflow origin a less likely scenario. Quantifying the velocity jump in Fig. 14 from the two peak positions in the position velocity diagram at 3.0 /-57.1 km s −1 and 6.1 /-49.2 km s −1 , we derive an estimate for the velocity gradient of ∼122 km s −1 pc −1 , more than twice the value found in the DCO + emission. These velocity jumps are discussed in Sect. 5.2 in the context of potential converging gas flows.
The temperature structure in ISOSS23053 shows a different distribution compared to what we saw before for ISOSS22478. Figure 15 shows the gas temperature derived again from the H 2 CO emission lines in comparison with the 1.3 mm continuum (top panel) and the integrated SiO emission (bottom panel). While we still find significantly large areas with temperatures around 10 K (dark blue in Fig. 15), similarly large areas are already above 30 K (lighter blue). However, even more importantly, ISOSS23053 exhibits, in significant parts of the map, temperatures above 50 K and reaching even 100 K and above. If we compare these areas of enhanced temperatures with the dust continuum emission, that traces the dense cores, and with the SiO emission, that should trace the shocked gas, we clearly find that the high-temperature gas is preferentially associated with the shocks. This is not only the case for the areas between the main cores, but we also find this temperature increases toward the stream-like SiO structures that likely trace outflowing gas. H. Beuther et al.: Kinematics of cloud collapse during high-mass star formation  However, we stress that the temperature enhancements as well as shock-tracing SiO emission are both found toward the region of the strong velocity jump in this regions. All these physical characteristics appear to be closely linked with what may happen in the case of converging gas flows (Sect. 5.2).

Fragmentation
Comparing the estimated masses and column densities to the corresponding values of the original CORE study of 20 highmass star-forming regions in the more evolved evolutionary stages of high-mass protostellar objects (HMPOs), massive young stellar objects (mYSOs), and ultracompact Hii regions (Beuther et al. 2018), we find that the masses are covering a similar range. In contrast to that, the column densities of the original CORE sources are between ∼10 23 and ∼10 25 cm −2 , about one order of magnitude larger than what we find here for the two younger regions. While this may be attributed to the younger evolutionary stage of the ISOSS sources presented here, we caution that the spatial resolution of the original CORE data was about a factor of 2 better (because of uniform weighting there and natural weighting here, see Sect. 3). Hence, with the higher spatial resolution, it was easier to resolve the highest column density peaks for the more evolved sources from the original CORE sample.
The mass-size relation shown in Fig. 4 is a bit puzzling because of the a priori uncertainty in the temperature estimates. While the T 20K approach results in a mass-size relation with M ∝ r 3.0±0.2 , which would correspond to constant mean volume densities for the sampled cores around 10 6 cm −3 , the T H2CO approach results in a flatter mass-size relation M ∝ r 2.0±0.3 . The latter agrees closely with the slope expected from Larson's 3rd relation (e.g., Ballesteros-Paredes et al. 2020) and corresponds to constant mean column densities. In the literature, mass-size relations with various exponents can be found. While some studies report relations close to the classical r −2.0 (e.g., Heyer et al. 2009;Lombardi et al. 2010;Ballesteros-Paredes et al. 2019), steeper relations have been found as well. For example, Kainulainen et al. (2011) found a mass-size relation with exponent 2.7 when sampling only the highest column density parts of their studied molecular clouds (those regions where the column density probability density functions leave the lognormal shape and rather flatten out). It is interesting to note that also in the recent comparison between the extreme environments of the nuclear starburst of the extragalactic system NGC253 and the central molecular zone of our Milky Way mass-size relations consistent with M ∝ r 3 were found (Krieger et al. 2020). Ballesteros-Paredes et al. (2012 discuss in detail the notion that mass-size relations with power-law indices from below 2 up to 3 have been reported in the literature (e.g., Mookerjea et al. 2004;Lada et al. 2008;Román-Zúñiga et al. 2010;Kainulainen et al. 2011;Könyves et al. 2015;Zhang & Li 2017;Veltchev et al. 2018). They reiterate that mass-size relations derived from constant column density thresholds (which is approximately the case for the 5σ intensity thresholds used here in the core finding algorithm), if the filling factor of the densest portions of the cores is small, should necessarily exhibit constant mean column densities, and thus, a M ∝ r 2.0 relationship. However, if the small filling-factor hypothesis for the densest gas in the continuum data and core identification were not fulfilled (e.g., Fig. 3), the cores may exhibit sharp boundaries and the lower column density material would not contribute substantially to the mean column density. Such observational conditions would imply that the column density probability density func- tions (PDFs) do not fall as fast as necessary to produce a M ∝ r 2.0 slope (Ballesteros-Paredes et al. 2012). Setting this picture into context with our observations, we cannot claim with certainty which temperature assumption is better for all cores and hence which mass-size relation, corresponding either to constant mean column density or to constant mean volume density. Possibly, some intermediate slope may represent the regions best. Hence, while our data are consistent with the classical 3rd Larson relation of constant mean column densities, some steeper slope, that results from not well sampled column density PDFs of the cores, cannot be excluded.
We can also set the derived masses and the projected nearestneighbor separations into context within the classical Jeans analysis by estimating the Jeans masses and Jeans lengths for the typical parental gas clumps. Assuming mean volume densities of the original gas clump of ∼ 10 5 cm −3 at temperatures of ∼20 K (e.g., Hennemann et al. 2008), the corresponding Jeans length and Jeans mass (e.g., Stahler & Palla 2005) are ∼17500 au and ∼0.9 M , respectively. Both of these values are approximately in Fig. 16. Core masses against nearest-neighbor separations from the minimum spanning tree analysis. The full pentagons show the results estimated assuming constant temperatures (T 20K ) whereas the crosses present the data assuming the temperatures derived from the H 2 CO line data (T H2CO ). The full line corresponds to the Jeans lengths and Jeans masses calculated at 20 K for a density grid between 5 × 10 3 and 10 6 cm −3 . For comparison, the dashed line corresponds to the Jeans lengths and Jeans masses calculated at a fixed density of 5 × 10 5 cm −3 with temperatures between 10 and 200 K. The error-bars in the top-right correspond to a factor 2 uncertainty for masses and 1 kpc distance uncertainty for the nearest-neighbor separations (resulting in the different error-bar end markers for the separations). the regime of the majority of core masses and separations (Table  2 and Fig. 5). Estimating the mean nearest-neighbor separations, we get values of ∼17170 au and ∼19560 AU for ISOSS22478 and ISOSS23053, respectively. These values also agree roughly with the estimated Jeans length, similarly to what was found by Sanhueza et al. (2019).
To show this global trend in more detail for the individual cores, Figure 16 presents the core masses versus the nearest neighbor separations, where the latter are used as proxies for the Jeans lengths. Corresponding relations for the Jeans analysis are shown there as well. While there is indeed scatter within the plot, the scatter is rather uniformly distributed around the typical Jeans relations. The scatter is slightly smaller for the masses derived with the T H2CO approach. Following Sanhueza et al. (2019), we can correct the nearest neighbor separations for average projection effects with an average correction factor of π/2 ∼ 1.57. This comparably small factor does not significantly affect the general correspondence of the Jeans relations and the data in Fig. 16. This is again different to the original CORE sample of more evolved sources where the derived parameters scattered above the Jeans relations ( Fig. 13 in Beuther et al. 2018). In that context, Beuther et al. (2018) considered what would happen if in the Jeans analysis, we replaced the thermal sound speed by the turbulent velocity dispersion. The expected relations would shift to the top-right in the mass-size relation, not improving the congruence between Jeans analysis and data. Furthermore, as visible in Figs. 7 and 11, the observed FWHM of individual spectral components are pretty low, typically around or even below 1 km s −1 (thermal line width of DCO + at 20 K ∼ 0.17 km s −1 ).
Hence, turbulent contributions to the line width are less strong than would be expected in the turbulent core scenario of highmass star formation (McKee & Tan 2003). Regarding filamentary fragmentation, in their ALMA study of very young highmass star-forming regions Svoboda et al. (2019) estimated the cylindrical fragmentation scale to be roughly a factor of 3.5 of larger than the typical Jeans length. In that cylindrical framework, the Jeans regimes outlined in Fig. 16 (dashed and continuous lines) would shift by that factor to the right, almost outside the plot. Hence, cylindrical fragmentation seems not the most important process in our regions either.
These considerations further support the idea that thermal Jeans fragmentation can explain our data and that turbulence as well as cylindrical fragmentation are unlikely to play a significant role in the fragmentation of at least these two very young high-mass star-forming regions. Similar results have also been reported in, for example, Palau et al. (2014Palau et al. ( , 2015Palau et al. ( , 2018, Henshaw et al. (2017), Cyganowski et al. (2017), Klaassen et al. (2017) as well as in the recent ALMA studies of very young high-mass star-forming regions by Svoboda et al. (2019) and Sanhueza et al. (2019). However, there are also reports in the literature that support the notion that for some other sources, turbulent pressure may be more important (e.g., Wang et al. 2011Wang et al. , 2014Zhang et al. 2015;Sadaghiani et al. 2020).

Gas dynamics
While the original CORE program focused more on the fragmentation of the parental gas clumps and the inner massive disks (Beuther et al. 2018), the CORE-extension program here focuses on the dynamics of the gas clumps. It is interesting to see that the observations of the gas dynamics of the two target regions do differ in many ways (Sect. 4.2).
The intermediate-mass region ISOSS22478 exhibits a more filamentary structure with velocity gradients along and across the filaments as well as multiple velocity components at the locations of the cores. Similar kinematic signatures along filamentary low-and high-mass star-forming regions have been reported in, e.g., Fernández-López et al. (2014), Henshaw et al. (2014), Beuther et al. (2015b) or Hacar et al. (2017a). The temperature distribution of ISOSS22478 exhibits relatively uniform low temperatures between roughly 10 and 30 K. Furthermore, we do not find strong signatures of shocks, neither in SiO emission nor specific temperature enhancements.
In comparison, the more massive region ISOSS23053 does not show a particularly filamentary structure, however, a few of the cores are aligned along a narrow spatial strip where two different velocity components of two separate cloud structures overlap. Position-velocity slices across that strip clearly reveal two velocity components and a very steep gradient between the two components of roughly 50 km s −1 pc −1 measured in DCO + (3 − 2) and even ∼122 km s −1 pc −1 measured in SiO (5 − 4) . This velocity gradient is significantly larger than that reported by Hacar et al. (2017a) of 5-7 km s −1 pc −1 for the gravitational collapse of the OMC-1 molecular cloud. Furthermore, for ISOSS23053, the region of velocity overlap shows strong SiO emission as well as a large temperature increase from roughly 20 K in the environmental cloud up to 100 K and more in the velocity overlap region. The combination of a velocity jump, enhanced SiO emission as well as temperature increase are all indicative of shocks that may have been caused by converging or colliding gas flows.
These two regions may at first sight appear to represent two potentially different modes of cloud formation. However, the dif-ferent signatures in the data could occur based on similar largescale cloud collapse processes. Multiple velocity components along filamentary structures have regularly been observed in various star-forming regions of different masses (e.g., Hacar et al. 2013Hacar et al. , 2018Kirk et al. 2013;Fernández-López et al. 2014;Henshaw et al. 2014;Tackenberg et al. 2014;Beuther et al. 2015b;Dhabal et al. 2018), but there are also significant observational differences reported. For example, some regions of low and high mass reveal sub-filaments (or fibers) that may collide and form larger-scale filamentary structures during the formation process (e.g., Hacar et al. 2013Hacar et al. , 2017bHacar et al. , 2018Henshaw et al. 2014;Smith et al. 2014). Other observations have revealed velocity gradients across filaments (e.g., Palmeirim et al. 2013;Fernández-López et al. 2014;Beuther et al. 2015b;Shimajiri et al. 2019) that may be explained by gravity-induced accretion from a magnetized sheet-like molecular cloud (e.g., Chen & Ostriker 2015;Chen et al. 2020). However, we should keep in mind that in these clouds, often not just one signature exists but that several signatures occur within the same region. For example, studies of Serpens south as well as of IRDC 18223 reveal velocity gradients perpendicular to the clouds in some part of the filaments, but in other parts of the clouds these studies find clearly distinct multiple velocity components that may resemble more the subfilament signatures (Fernández-López et al. 2014;Dhabal et al. 2018;Beuther et al. 2015b). Chen et al. (2020) point out that multiple velocity-components along a filament do not necessarily always have to be real structures but that the observational bias that different spectral lines trace different density regimes can also create artificial multiple velocity components in observations (see also Ballesteros-Paredes & Mac Low 2002;Zamora-Avilés et al. 2017). Different inclination angles and optical depth effects may also induce a similar bias in some data as well.
From a physical point of view, both scenarios, namely, colliding pre-existing filaments and gravity-induced accretion in sheet-like structures have a feature in common: both basically rely on the large-scale gravitational collapse of the molecular clouds. Low magnetization may result in more early fragmentation and then interacting pre-existing sub-filaments whereas higher magnetization should result in smoother structures and more coherent accretion onto filaments (see for example, Hennebelle & Inutsuka 2019, and references there in).
The distinct signature of a steep velocity gradient or even sharp velocity jump in ISOSS23053 can be interpreted as a signature of converging and colliding gas flows. For example, the cloud formation models with colliding flow conditions by Gómez & Vázquez-Semadeni (2014) produce similar velocity gradient structures, as seen in ISOSS23053. However, other simulations of cloud collapse, either hydrodynamic or magnetohydrodynamic, can produce qualitatively similar velocity gradients as well (e.g., Smith et al. 2013;Chen et al. 2020). The outstanding aspect of the ISOSS23053 case is the steepness of the gradient with ≥50 km s −1 pc −1 . While variations in inclination angles can induce differences in the slope of such a gradient, it is also possible that different physical properties, for instance, different magnetization or external compression may enhance the observed gradients. It appears that the smaller and smoother velocity gradients, as observed, for example, in the Serpens south, IRDC 18223, or the Orion integral shape filament (Fernández-López et al. 2014;Beuther et al. 2015b;Hacar et al. 2017a), can be explained mainly by gravitational collapse, whereas the sharp velocity jump observed here in ISOSS23053 may require an external cause that drives the converging gas flow. Gravitational collapse can then occur at the colliding flow interface.
Differentiating these variations requires additional observational as well as theoretical work. For example, it will be important to infer the magnetic field for a sample of clouds and investigate whether steep velocity gradients or multiple velocity components may correlate with the strength or orientation of the environmental magnetic fields. One obvious difference between the two regions is their Galactic latitude. While ISOSS23053 is close to the Galactic plane with a latitude of −0.28 deg, ISOSS22478 is located at larger Galactic latitudes of +4.26 deg. It is reasonable that closer to the Galactic plane, along with more star formation activity, there is also a more dynamic environment that can further amplify the potential dynamic gas structure and velocity gradients. On the simulation side, we would need to sample a broader range of magnetizations to infer whether stronger magnetic fields indeed favor smoother accretion onto filamentary structures, while lower magnetization may favor early formation of sub-filaments that could then collide during ongoing collapse.
What all of the above scenarios have in common is that they are based on a dynamical cloud collapse scenario that forms filamentary structures during the collapse motions. Environmental effects such as external shocks, converging gas flows or different magnetizations come into play, but they are all part of a dynamical collapse and fragmentation scenario.

Conclusions and summary
In the context of the NOEMA large program CORE, which studies, among other topics, the fragmentation processes during high-mass star formation, we studied two young regions of intermediate-to high-mass star formation, namely ISOSS22487 and ISOSS23053. We used NOEMA and the IRAM 30 m telescope to produce high spatial resolution (∼0.8 ) mosaics to investigate their larger-scale fragmentation and cloud formation processes. In both regions, many fragments can be identified in the 1.3 mm continuum emission, and a diversity of spectral lines can be imaged over the entire field of view. We concentrate on the fragmentation based on the continuum emission, and on the kinematics of the clouds based on the DCO + , H 2 CO, and SiO emission. A forthcoming paper will present the entire spectral line data and a chemical analysis (Gieser et al. in prep.). The 1.3 mm continuum data reveal 15 and 14 cores in ISOSS22478 and ISOSS23053, respectively. Most of them are associated with filamentary sub-structures in the cloud. Depending on the assumed temperatures of the cores, we estimate slightly different core masses. These differences then result in mass-size relations varying possibly between M ∝ r 2 and M ∝ r 3 . Hence, while the data are consistent with the classical mass-size relation of r 2 expected from the 3rd Larson relation for regions with constant mean column densities, a steeper relation cannot be ruled out. The latter could be explained because the cores in these maps have sharp boundaries and, thus, the lower iso-contours in the mm emission do not have much larger areas than higher iso-contours. Such steeper mass-size relation would correspond to similar mean densities for the cores. Furthermore, the correlation of the core masses with their nearest-neighbor separations is consistent with thermal Jeans fragmentation. In contrast to the original CORE sample of slightly more evolved sources, here, we barely find separations at our spatial resolution limit. Hence, our data resolve the large-scale fragmentation of the parental gas cloud well.
The kinematic analysis reveals very different observational signatures between the two regions. While ISOSS22478 is more filamentary in nature with several velocity components along the length of the filament, ISOSS23053 shows at the locations of the most massive cores a very steep velocity jump between two velocity components. The velocity gradient is roughly 50 km s −1 pc −1 measured in DCO + (3 − 2) and even ∼122 km s −1 pc −1 measured in SiO (5 − 4). While these signatures appear disjointed at first sight, we argue that all these kinematic observations can be understood within the framework of dynamically collapsing clouds. Depending on the environmental properties, for example, external shocks causing converging gas flows or magnetized sheets that preferentially form filaments that accrete gas from the parental gas structure, different observational signatures of a global cloud collapse may be present. Further observational and theoretical investigations are needed to pin down in more detail how much the physical processes (e.g., magnetization, external compression) as well as the observational biases (e.g., tracers of different density, optical depth effects, inclination angles) contribute to the overall dynamical collapse of the clouds.