© C. Gieser et al. 2021

1 Introduction

Massive stars (M_⋆ ≥ 8M_⊙) are rare, compared to low-mass stars with M_⋆ ≲ 1M_⊙, but important in shaping the physical and chemical properties of a galaxy. For example, their high bolometric luminosities dominate the stellar contribution to the spectral energy distribution (SED) of galaxies (e.g., Walcher et al. 2011) and thus provide a high flux of ionizing radiation. Energetic, mechanical, and radiative feedback is provided to the interstellar medium (ISM) through outflows (e.g., Beuther et al. 2002b; Kölligan & Kuiper 2018), stellar winds (e.g., Meynet et al. 1994; Gatto et al. 2017), and supernova explosions (e.g., McKee & Ostriker 1977; Girichidis et al. 2016).

In the earliest stages of high-mass star formation (HMSF), the stellar birthplaces are still deeply embedded within their dense parental molecular cloud and not visible at optical or infrared (IR) wavelengths. Therefore, the cold gas and dust can only be studied at millimeter (mm) or sub-mm wavelengths. For comprehensive reviews of massive star formation, we refer to, for example, Beuther et al. (2007a), Bonnell (2007), Zinnecker & Yorke (2007), Smith et al. (2009), Tan et al. (2014), Schilke (2015), Motte et al. (2018), and Rosen et al. (2020).

Once a protostar forms, gas accretion and bipolar outflows are commonly observed, and the young stellar objects (YSOs) become bright at mid-infrared (MIR) wavelengths. HMSF is generally clustered and observed toward giant molecular clouds with sizes of R ≳ 1 pc (e.g., Roman-Duval et al. 2010), providing a large mass reservoir for individual massive clumps with typical sizes of R ≳ 0.1 pc. Following the nomenclature of Beuther et al. (2007a) and Zhang et al. (2009), these clumps are commonly observed to fragment into individual cores at high angular resolution, R ≈ 0.01−0.1 pc, where a single or a small system of multiple protostars forms.

A large diversity of core fragmentation properties is observed. In some regions, there is a single dominant massive core (e.g., Beuther et al. 2018; Maud et al. 2019), while in other regions a large number of cores are found (e.g., Palau et al. 2014, 2015; Beuther et al. 2018, 2021). In addition, it has been suggested that HMSF could occur at hubs of filamentary structures (e.g., Myers 2009; Galván-Madrid et al. 2010, 2013; Schneider et al. 2012; Tigé et al. 2017; Kumar et al. 2020). It is not yet fully understood what drives the accretion of the surrounding gas and dust from cloud and clump scales onto the cores (e.g., Smith et al. 2009). Variations in magnetic field strengths (e.g., Hennebelle & Inutsuka 2019) and density profiles could also be important factors that influence the multiplicity and core mass (e.g., Commerçon et al. 2011; Girichidis et al. 2011).

High-mass star-forming regions (HMSFRs) are sites with a rich molecular content and for dedicated reviews we refer to Herbst & van Dishoeck (2009) and Jørgensen et al. (2020). Most of the total of ~200 molecules detected in the ISM are found toward HMSFRs (McGuire 2018). Chemical reactions occur in the gas phase, on dust grain surfaces, and in the icy mantle layers (e.g., Garrod et al. 2008). A current challenge is to understand, for example, the chemical segregation between nitrogen(N)- and oxygen(O)-bearing species (e.g., Wyrowski et al. 1999; Jiménez-Serra et al. 2012; Feng et al. 2015; Allen et al. 2017; Gieser et al. 2019).

In order to study HMSF at core scales, we carried out the survey “CORE - Fragmentation and disk formation during high-mass star formation”, which is a large program with the NOrthern Extended Millimeter Array (NOEMA) targeting 18HMSFRs in the northern hemisphere with the currently highest possible angular resolution at 1.3 mm (~0.″4). An overview of the CORE project and analysis of the 1.3 mm continuum data are presented in Beuther et al. (2018). The molecular line data, revealing a high degree of chemical complexity on scales <10 000 au, are analyzed by Gieser et al. (2021).

The observations of the CORE sample were carried out between 217.2 and 220.8 GHz (~1.3 mm) with the WideX correlator and a single pointing toward each region. Two additional pilot regions were previously observed in a similar setup (Feng et al. 2016).

In the physical and chemical analysis of the full CORE sample (Gieser et al. 2021), 22 objects are analyzed, defined as “cores”, that have a radially decreasing temperature profile and associated 1.3 mm continuum emission. Molecular column densities are derived from the 1.3 mm spectra (e.g., C¹⁸O, SO, OCS, HNCO, HC₃N, CH₃OH, CH₃CN), and each core is modeled with the physical-chemical model MUlti Stage ChemicaL codE (MUSCLE) to estimate chemical timescales, τ_chem. The luminosities toward these regions are high (L > 10⁴ L_⊙; Beuther et al. 2018), suggesting that these regions contain massive YSOs and that the evolutionary stages of most of the regions are between the high-mass protostellar object (HMPO) and hot molecular core (HMC) phase. Some regions have already formed hyper- or ultra-compact HII (HCHII or UCHII) regions. For a more detailed description of the physical properties of the IRDC, HMPO, HMC, and UCHII evolutionary phases, we refer to, for example, Gerner et al. (2014) and Gieser et al. (2021). With MUSCLE, gradients of τ_chem are found for cores located within the same region if the separations are large, while nearby cores have similar chemical timescales.

In this paper, we study the core properties in the two young and cold intermediate- to high-mass star-forming regions ISOSS J22478+6357 and ISOSS J23053+5953. The data were obtained as an extension of the CORE project (“CORE-extension”) with both regions observed in multiple pointings creating large mosaics and the new PolyFiX correlator increasing the spectral bandwidth by a factor of four. An overview and the analysis of the fragmentation and kinematic properties is presented in Beuther et al. (2021). These two additional regions are at very early evolutionary stages, similar to typical infrared dark clouds (IRDCs), with T_clump ≈ 20 K (Ragan et al. 2012). Thus, these two sources complement the CORE sample with regions that tend to occupy an earlier phase in the evolutionary sequence of HMSF (Beuther et al. 2007a). Both regions were selected based on targets of the ISOPHOT Serendipity Survey (ISOSS) observing the sky during slew times at 170 μm (Krause et al. 2004).

The ISOSS J22478+6357 region was first identified with the Infrared Astronomical Satellite (IRAS) as IRAS 22460+6341 (Beichman et al. 1988), and is also known as J224749.9+635647 (Di Francesco et al. 2008). We note that the associated IRAS source (Kerton & Brunt 2003) is spatially offset by ~41″ from the studied IR and mm counterparts in this and other more recent studies (Hennemann et al. 2008; Ragan et al. 2012; Bihr et al. 2015; Beuther et al. 2021). Wide-field Infrared Survey Explorer (WISE) images at 12 and 22 μm show extended emission around the IRAS source, which can be interpreted as a signpost of intermediate-mass star formation (Lundquist et al. 2014).

With Herschel observations using the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) instrument as part of the Earliest Phases of Star Formation (EPoS) program, Ragan et al. (2012) determine a kinematic distance of $d={3.23}_{-0.60}^{+0.61}$ kpc (using the model by Reid et al. 2009) and a total gas mass of 104 M_⊙ within a total of seven resolved clumps. Thesystemic velocity v_LSR is − 39.7 km s⁻¹. The region is located in the outer galaxy, d_gal = 9.7 kpc away from the Galactic center.

A detailed study of a few ISOSS regions, including ISOSS J22478+6357, is presented in Hennemann et al. (2008). These authors detect two submm clumps with the Submillimetre Common User Bolometer Array (SCUBA) at the James Clerk Maxwell Telescope (JCMT) at 450 and 850 μm, “SMM1 E” and “SMM1 W”, with sizes of 0.14 and 0.24 pc, gas masses of 64 M_⊙ and 116 M_⊙, and dust temperatures of 15 and 14 K, respectively. By fitting the SED of a bright MIR source, one of the more evolved YSOs in the SMM1 E clump is classified as an intermediate-mass star with central mass M_⋆ = 6−8.5M_⊙ and a system age of 10⁶−6 × 10⁶ yr. Using the 850 μm JCMT SCUBA observations, Bihr et al. (2015) estimate a total gas mass of 140 M_⊙ in ISOSS J22478+6357.

The region ISOSS J23053+5953 was originally identified as a point source by IRAS under the identifier IRAS 23032+5937 (Beichman et al. 1988). In the literature it is also known as J230523.6+595356 (Di Francesco et al. 2008) or G109.995-00.282 (Rosolowsky et al. 2010). With the Bolocam Galactic Plane Survey (BGPS), a 1.1 mm flux of 0.826 Jy is derived within a radius of R = 38.″32 (Rosolowsky et al. 2010). Schlingman et al. (2011) estimate the following properties using the BGPS data: M = 130M_⊙, beam average density n = 3000 cm⁻³, volume averaged density n = 840 cm⁻³, free-fall timescale τ_ff = 1.4 × 10⁶ yr, and crossing timescale τ_cross = 4.0 × 10⁵ yr.

Single-dish observations with a linear resolution of ≈0.3−1 pc reveal molecular line emission toward the region with line widths of 1.5− 5 km s⁻¹ (e.g., CO, CS, HCO⁺, N₂ H⁺, CH₃ CCH, CH₃ CHO; Harju et al. 1993; Wouterloot et al. 1993; Bronfman et al. 1996; Alakoz et al. 2002; Shirley et al. 2013; Vasyunina et al. 2014). Based on NH₃ observations with the Effelsberg 100 m telescope, the kinetic temperature was estimated to be ~15 K (Wouterloot et al. 1988; Harju et al. 1993). Vasyunina et al. (2014) estimate a kinetic temperature of 31 ± 1 K based on CH₃ CCH line emission. Early studies found that there are two velocity components in molecular line emission, for example, for NH₃ (at − 52.3 and − 50.8 km s⁻¹, Wouterloot et al. 1988), H₂O (at − 52.3 and − 49.5 km s⁻¹, Wouterloot et al. 1993), and CS (at − 51.8 and − 51.2 km s⁻¹, Larionov et al. 1999).

Very Large Array (VLA) and Effelsberg 100 m observations of the NH₃ (1,1) and (2,2) lines at an angular resolution of 4″ are analyzed by Bihr et al. (2015). These authors find a steep velocity gradient of >30 km s⁻¹ pc⁻¹ toward the region suggesting a dynamical collapse and/or converging gas flow. With JCMT SCUBA 850 μm observations, a total gas mass of 610 M_⊙ is estimated by the authors in ISOSS J23053+5953. The kinematic data of the CORE-extension project further resolve this steep velocity gradient in DCO⁺ (3− 2) being higher than 50 km s⁻¹ pc⁻¹ (Beuther et al. 2021).

Wouterloot & Walmsley (1986) and Wouterloot et al. (1993) detected H₂O maser emission. But no H₂O maser emission was detected in follow-up studies (Comoretto et al. 1990; Palagi et al. 1993; Palla et al. 1993; Slysh et al. 1999; Valdettaro et al. 2001; Sunada et al. 2007), so the presence or potential variability of H₂O maser emission remains unclear. No CH₃OH maser emission was detected (Wouterloot et al. 1993; Kalenskii & Sobolev 1994).

ISOSS J23053+5953 is a target of the EPoS survey as well (Ragan et al. 2012). The authors derive a gas mass of 488 M_⊙ within three resolved clumps and a kinematic distance of $d={4.31}_{-0.62}^{+0.64}$ kpc at a systemic velocity of − 52.5 km s⁻¹ (using the model by Reid et al. 2009). This is in agreement with other distance estimates, for example, by Yang et al. (2002, d = 4.82 kpc), Schlingman et al. (2011, d = 4.253 kpc), Ellsworth-Bowers et al. (2015, $d={4.44}_{-0.70}^{+0.74}$ kpc), and the Planck Collaboration XXVIIII (2016, d = 5.07 kpc). With a galactocentric distance of d_gal ≈ 10 kpc (Wouterloot & Brand 1989; Ellsworth-Bowers et al. 2015) the region is located in the outer galaxy as well.

Birkmann et al. (2007) find a young and accreting massive protostar toward one of the clumps with an associated outflow. The line profile of the optically thick HCO⁺ 3− 2 transition indicates infalling material with self-absorption in the red-shifted line wing compared the optically thin H¹³CO⁺ 3− 2 isotopologue.

Pitann et al. (2011) studied the region using Spitzer observations with the InfraRed Array Camera (IRAC, Fazio et al. 2004) and the Multiband Imaging Photometer for Spitzer (MIPS, Rieke et al. 2004) instruments. These authors find that the region, harboring a cluster of sources, has extended polycyclic aromatic hydrocarbons (PAH) emission and warm dust components. In addition, the emission of forbidden transitions of FeII, SiII, and SI imply post-shocked gas that occurred by a J-shock. This is confirmed by the detection of high-energy H₂ transitions from S(0) to S(7) with excitation energies ranging between 510−7200 K. These authors do not find an indication of a photodissociation region (PDR).

An overview of the 1.3 mm observations with NOEMA, the core fragmentation properties, and a detailed analysis of the kinematic properties of ISOSS J22478+6357 and ISOSS J23053+5953 are presented in Beuther et al. (2021). We expand the analysis with a detailed investigation of the molecular lines detected at 1.3 mm. In addition, we used archival MIR and far-infrared (FIR) data to study the clustered nature of both regions. The MIR observations reveal emission of more evolved protostars while the FIR emission traces the cold dust emission.

The paper is organized as follows: The NOEMA and Institut de RadioAstronomie Millimétrique (IRAM) 30 m telescope observations and data calibration are described in Sect. 2. The analysis of the continuum data is given in Sect. 3. The molecular line data are analyzed in Sect. 4. In Sect. 5 we apply a physical-chemical model to a few line-rich cores. We discuss our results in Sect. 6 and a summary and conclusions are given in Sect. 7.

2 Observations

2.1 CORE-extension data

The NOEMA mosaic observations were taken in February and March 2019 with ten antennas in the A, C, and D array configurations at 1.3 mm (Band 3) using the PolyFiX correlator. ISOSS J22478+6357 and ISOSS J23053+5953 were observed with six and four NOEMA pointings, respectively. A summary of both regions, for example, the coordinates and velocities, is shown in Table 1. Complementary IRAM 30 m observations in the same frequency range, to include short-spacing information of the spectral line data, were obtained using the Eight MIxer Receiver (EMIR, Carter et al. 2012) in June 2019.

2.1.1 NOEMA data

The NOEMA data were calibrated using the CLIC package in GILDAS¹. The PolyFiX correlator simultaneously covers ~8 GHz in two sidebands (lower sideband, LSB, and upper sideband, USB) and in the two orthogonal linear polarizations (horizontal and vertical) with a fixed channel spacing of 2 MHz (~2.7 km s⁻¹ at 1.3 mm). Rest frequencies from 213.3 GHz to 221.3 GHz in the LSB and from 228.7 to 236.7 GHz in the USB were covered by the observations. High-resolution units with a channel spacing of 62.5 kHz (~0.084 km s⁻¹ at 1.3 mm) were placed within the broadband correlator units.

The NOEMA 1 mm data of the original CORE sample consisted of single pointing observations toward each region, which allowed us to self-calibrate these data sets with the GILDAS selfcal task (a detailed description of the method is presented in Gieser et al. 2021). Since self-calibration of mosaic observations is currently not possible, we use the standard calibrated NOEMA data for the CORE-extension data presented in Beuther et al. (2021) and in this work.

2.1.2 IRAM 30 m data

The EMIR data were calibrated using the CLASS package in GILDAS. The EMIR instrument consists of four basebands with a width of ~4 GHz in each baseband: the lower outer (LO, 213.5−217.5 GHz), lower inner (LI, 217.3−221.3 GHz), upper inner (UI, 229.2−233.2 GHz), and upper outer (UO, 233.0−237.0 GHz) baseband. The chosen Fast Fourier Transform Spectrometer (FTS) backend delivers a channel separation of 200 kHz (~0.27 km s⁻¹ at 1.3 mm). The half power beam width (HPBW) is 11.″8 in the lower basebands and 11.″2 in the upper basebands. After an initial inspection of the data, spiky channels caused by noise artifacts were filled with Gaussian noise. The antenna temperature $T_{\text{A}}^{*}$ was converted to main beam temperature T_mb using $T_{\text{mb}}=\frac{F_{\text{eff}}}{B_{\text{eff}}}\times T_{\text{A}}^{*}$ with F_eff = 0.915 and B_eff = 0.57.

2.1.3 Imaging

The NOEMA continuum and combined NOEMA + IRAM 30 m (“merged”) spectral line data were imaged using the MAPPING package in GILDAS. Primary beam correction was applied to the final continuum and spectral line data products.

In a first inspection of the NOEMA spectral line data of both regions, we identified all detected emission lines. An overview of all lines and their properties is shown in Table 2. All of these emission lines were covered by the IRAM 30 m EMIR observations, and hence the interferometric and single-dish data can be combined in order to recover missing flux filtered out by the interferometer. If available, we used the data obtained with the high-resolution units smoothed to a spectral resolution of 0.5 km s⁻¹ to increase the signal-to-noise ratio (S/N) and otherwise with the low-resolution basebands smoothed to a spectral resolution of 3.0 km s⁻¹ to increase the S/N.

The continuum is extracted from the low-resolution spectral line data by masking out all channels with line emission using the uv_filter and uv_continuum tasks. The continuum data in the LSB and USB were merged with the task uv_merge. The continuum data were CLEANed with natural weighting using the Clark algorithm (Clark 1980). The synthesized beam (θ_maj × θ_min, PA) and noise σ_cont of the continuum image are 0.″92 × 0.″73, 50° and 0.057 mJy beam⁻¹ for ISOSS J22478+6357 and 0.″84 × 0.″74, 67° and 0.16 mJy beam⁻¹ for ISOSS J23053+5953. Short-spacing information is only available for the spectral line data, so spatial filtering affects the continuum data. In Beuther et al. (2021) it is estimated that only 10− 20% of the total flux is recovered by the NOEMA observations.

The continuum was subtracted from the spectral line data with the task uv_baseline. For each visibility, a baseline is fitted and subtracted with all channels of line emission masked out (listed in Table 3). The NOEMA observations were combined with the IRAM 30 m data using the uvshort task.

The continuum-subtracted spectral line data were CLEANed with natural weighting using the SDI algorithm (Steer et al. 1984) for lines with extended emission within the field-of-view (all CO isotopologues, SO, SiO, and H₂CO) and with the Clark algorithm (Clark 1980) for lines with less extended emission. A detailed comparison of these two CLEAN algorithms applied to the CORE region W3 IRS4 is presented in Mottram et al. (2020). The properties of the spectral line data products are summarized in Table 3. The line noise σ_line, computed in emission-free channels, is ~0.2 and ~0.07 K channel⁻¹ in the high-resolution and low-resolution line data, respectively.

2.2 Archival data

In addition to the high angular resolution NOEMA data at 1.3 mm, we used archival MIR to FIR observations obtained with the Spitzer and Herschel space telescopes to study the clustered nature of the regions and constrain the SEDs. For both regions, archival data covering the full field-of-view (FOV) of the NOEMA mosaic exists.

In the FIR, both regions were targets in the Herschel key program EPoS and the properties of the cold dust are analyzed in Ragan et al. (2012). Three photometric bands of the PACS instrument at 70, 100, and 160 μm, with angular resolutions of ~5″, ~7″, and ~12″, respectively, cover the peak of the SED of cold dust emission. Science-ready data products were taken from the Herschel Science Archive.

The Spitzer IRAC instrument has four photometric bands at 3.6, 4.5, 5.8, and 8.0 μm. The angular resolution is 1″−2″, so comparable to the 1.3 mm observations with NOEMA. The Spitzer MIPS observations at 24 μm have an angular resolution of 7″. The 70 μm MIPS data were not used, as the Herschel PACS 70 μm data have a higher angular resolution. Science-ready data products were obtained from the Spitzer Heritage Archive. For a detailed description of the Spitzer observations, we refer to Hennemann et al. (2008) for ISOSS J22478+6357 and to Birkmann et al. (2007) and Pitann et al. (2011) for ISOSS J23053+5953.

3 Continuum

The NOEMA 1.3 mm continuum data reveal the compact dust emission in the two star-forming regions. The subarcsecond resolution achieved with NOEMA allows us to study individual fragmented millimeter (mm) cores. In this section, the core properties are analyzed using the NOEMA 1.3 mm continuum emission in combination with archival MIR and FIR data.

3.1 Fragmentation properties

The NOEMA 1.3 mm continuum data are shown in Fig. 1 for ISOSS J22478+6357 and ISOSS J23053+5953. In both regions, the continuum data reveal a large number of fragmented mm cores. In total, ten bright mm cores with S∕N > 20 are detected, four in ISOSS J22478+6357 and six in ISOSS J23053+5953.

The fragmentation properties are studied in detail in Beuther et al. (2021). By applying the clumpfind algorithm (Williams et al. 1994)to the continuum data, 15 and 14 individual mm cores can be identified for ISOSS J22478+6357 and ISOSS J23053+5953, respectively.The cores are labeled in Fig. 1. In this study, we investigate the physical and chemical properties of these 29 cores. Anoverview and summary of the results from Beuther et al. (2021), including the core position, 1.3 mm peak intensity I_peak, 1.3 mm integrated flux F₁₃₀₀, and outer radius r_out derived with clumpfind, is shown in Table 4.

A comparison of the NOEMA 1.3 mm continuum and archival MIR and FIR observations obtained with the Spitzer and Herschel space telescopes is presented in Fig. 2. The Spitzer IRAC and MIPS observations highlight evolved protostars. The Herschel PACS observations at 70, 100, and 160 μm trace the cold dust emission on clump scales. The mm cores have counterparts not only at FIR wavelengths, such that they are embedded in the cold clumps, but there are also a few cores that are MIR-bright (e.g., mm core 13 in ISOSS J22478+6357 and mm core 2 in ISOSS J23053+5953). The SED of mm cores with associated Spitzer IRAC sources are constructed and modeled in Sect. 3.3.

Based on the 1.3 mm continuum peak intensity I_peak and integrated flux F₁₃₀₀ of each core, Beuther et al. (2021) calculate the molecular hydrogen column density N(H₂) and core mass M_core assuming optically thin dust emission. For the kinetic temperature T_kin, the H₂CO rotation temperature T(H₂CO) was used assuming local thermal equilibrium (LTE) conditions (Sect. 4.4). In order to check if the assumption of optically thin dust emission is valid, we computed the 1.3 mm continuum optical depth τ_cont at the position of the peak intensity of each core with $${\tau }_{\text{cont}}=-\text{ln}(1-\frac{I_{\text{peak}}}{B_{\nu }(T_{\text{kin}})}),$$(1)

where B_ν is the Planck function. The continuum optical depth, τ_cont, as well as N(H₂), M_core, and T_kin for each core are summarized in Table 4. The optical depth is <0.01 for all 29 cores; therefore, the assumption of optically thin dust emission is valid.

The core masses M_core, derived from the integrated 1.3 mm flux, vary between 0.04 and 4.61 M_⊙ (Table 4). This indicates that most of the mm cores form low- and intermediate-mass stars. The sum of all core masses is 17 M_⊙ and 23 M_⊙ in ISOSS J22478+6357 and ISOSS J23053+5953, respectively. The mass estimates using single-dish observations (Ragan et al. 2012; Bihr et al. 2015) are a factor of 10−30 higher than estimated from the 1.3 mm NOEMA data. The interferometric observations of the dust continuum suffer from spatially filtering the extended emission, the estimated core masses are therefore only lower limits. The single-dish observations at FIR wavelengths from Herschel PACS and JCMT SCUBA reveal that both regions have a large gas mass reservoir of a few 100 M_⊙ (Ragan et al. 2012; Bihr et al. 2015). This provides a mass reservoir for further growth and indicates that mass assembly is not complete. Bright MIR sources suggest the presence of at least intermediate-mass protostars in both regions.

Birkmann et al. (2007) derive core masses of 26, 4.4, and 4.4 M_⊙ for cores 1, 2, and 3, respectively, in ISOSS J23053+5953 using 1.3 mm observations. These values are higher than the estimates by Beuther et al. (2021). As in the calculation by Birkmann et al. (2007) the assumed kinetic temperature was derived from lower dust temperatures based on the single-dish observations. The estimated core mass of mm core 13 in ISOSS J22478+6357 is low (0.44 M_⊙), but Hennemann et al. (2008) infers that this MIR source (their “source 3”) has a central mass of 6−8.5M_⊙ using the YSO SED models by Robitaille et al. (2007). By applying the same SED models to mm core 2 (their “SMM East”) in ISOSS J23053+5953, Pitann et al. (2011) estimate a central mass of 4.1−6.6M_⊙, whereas Beuther et al. (2021) derive a core mass of 2.51 M_⊙.

The calculation of the core mass depends strongly on the assumed temperature. In Beuther et al. (2021), core masses are also derived assuming a uniform temperature of 20 K. In cases when the kinetic temperature is actually higher, derived from the H₂CO rotation temperature maps (Sect. 4.4), the core masses are overestimated by a factor of a few. While the assumption of a low temperature at clump scales is valid, the high angular resolution data reveals that some mm cores already harbor protostars that heat up the surrounding envelope. Differences in the mass estimates are also caused by a higher assumed gas-to-dust mass ratio (150, Beuther et al. 2021) instead of the canonical value of 100 used in the literature. Given the mass reservoir and that the regions are in an early evolutionary stage and thus gas accretion will continue for a few 10⁵ −10⁶ yr, some of the low-mass cores might form intermediate- to high-mass protostars in the future.

Fig. 1 NOEMA 1.3 mm continuum emission of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The NOEMA continuum data are presented in color scale and pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. The core outer radius rout is indicated by a dashed blue circle. The synthesized beam is shown in the bottom left corner. A linear spatial scale of 0.1 pc is indicated by a black scale bar.

3.2 Radial density profiles

Observations of individual high-mass cores suggest that the density profiles of the gas and dust envelope around the forming protostars can be described by power-law profiles with power-law indices p typically ranging from 1.5 to 2.0 (e.g., Beuther et al. 2007b; Zhang et al. 2009; Palau et al. 2014, 2021; Gieser et al. 2021; Gómez et al. 2021). The density profiles are described as $$n(r)=n_{\text{in}}\times {(\frac{r}{r_{\text{in}}})}^{-p},$$(2)

with density n_in and arbitrary radius r_in. The density profiles can be inferred, for example, from the radial intensity profiles of the dust continuum emission (van der Tak et al. 2000; Beuther et al. 2002a; Palau et al. 2014). While the aforementioned studies are based on single-dish observations, the interferometric continuum data suffer from significantly spatially filtering the extended emission (Beuther et al. 2021). This effect can cause intensity profiles in the image domain to be much steeper than in reality and result in too steep (p > 3) density profiles. In order to overcome this problem, we therefore analyze the continuum data in the uv plane - excluding the regime of missing baselines - in order to estimate the power-law index p of each core. The same method to infer the core density profiles as presented in Gieser et al. (2021) is applied to the mm cores. Assuming optically thin dust emission (which is valid for all cores, Table 4), the power-law index p can be calculated from the power-law index of the azimuthally averaged visibility amplitudes $\overline{V}~s^{\alpha }$ (taking into account both the real and imaginary components), with uv distance s (e.g., Adams 1991; Looney et al. 2003). The density power-law profile then equals to $$p=\alpha -q+3,$$(3)

where q is the power-law index of the temperature profile. The temperature profile T(r) of the cores is also assumed to have a spherically symmetric distribution and is described as $$T(r)=T_{\text{in}}\times {\left(\frac{r}{r_{\text{in}}}\right)}^{-q},$$(4)

with temperature T_in = T(r_in) and arbitrary radius r_in. The temperature profiles are estimated from the H₂CO temperature maps in Sect. 4.5.

Many of the fainter mm cores are only marginally or not resolved. In order to reliably derive the radial density and temperature profiles, we restrict the analysis to mm cores that are detected with S∕N > 20 in the 1.3 mm continuum data. This is the case for cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953.

The azimuthally averaged visibility amplitude $\overline{V}$ and its standard deviation were computed for each core with the uvamp task in MIRIAD (Sault et al. 1995) considering the real and imaginary components.The phase center was shifted to the location of the core (Table 4). The visibilities were binned in a bin size of 20 up to 600 kλ, which corresponds to the longest NOEMA baseline (774 m).

The computation of the azimuthally averaged visibility amplitudes of a certain core can be influenced by other mm-bright cores in the field (e.g., in ISOSS J23053+5953 mm core 2 is close to mm core 3, Fig. 1). Thus, before the radial visibility amplitudes of a certain core were computed, we first subtracted the emission of the remaining four brightest cores in the region. Using the uv_fit task in GILDAS MAPPING, we found that the emission of the mm cores is modeled best by fitting two circular Gaussian functions with a full width half maximum (FWHM) of 0.″5, and 2.″0, respectively,in order to take into account both compact and extended emission. The flux of the two components was varied for each core separately in order to take into account that the cores have a varying morphology (e.g., compact core or a compact core embedded in a more extended envelope). The model emission of the four remaining brightest mm cores, described each by the two Gaussian functions, was then subtracted from the original uv data set. We carefully checked, for each core, that the emission of the four brightest remaining cores was subtracted correctly by imaging the cores-subtracted uv data set.

As an example, after core 3 in ISOSS J23053+5953, the four remaining brightest cores within the FOV are cores 1, 2, 4, and 5 (Table 4). We first fitted and subtracted the model emission of cores 1, 2, 4, and 5 with the uv_fit task. We imaged this cores-subtracted uv data set and validated that the emission of cores 1, 2, 4, and 5 was removed without affecting the emission of core 3. Then we shifted the phase center in the cores-subtracted uv data to the location of core 3 and computed the azimuthally averaged visibility amplitudes using the uvamp task as described above.

The azimuthally averaged visibility profiles for cores 1− 4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953 are shown in Fig. 3. The radial visibility profiles of most of the cores follow a single power-law profile. Only at the shortest baselines, at 0− 20 kλ, the profiles flatten for most of the resolved cores. The flattening at short uv distances is caused by spatial filtering since the shortest baseline is ~18 m. In general, at large uv distances (> 300 kλ), there is an increased scatter of the data points due to the fact that the number of long baselines of the NOEMA interferometer is smaller compared to the number of short baselines. This also causes an increase in the uncertainties of the binned data points. Cores 3 and 4 in ISOSS J22478+6357 have a flat profile, since these cores are only marginally resolved (Fig. 1). The visibility profile of core 1 in ISOSS J23053+5953 flattens at uv distances <60 kλ. A large bipolar molecular outflow is observed toward this core (Sect. 4.1) that might impact the radial profile of the dust envelope.

Cores 3 and 4 in ISOSS J22478+6357 and cores 2 and 3 in ISOSS J23053+5953 are very close, 1.″2 (3900 au) and 2.″9 (12 000 au), respectively (Fig. 1). The subtraction of the emission of the nearby cores is difficult in these cases, since the cores are embedded within a common envelope.

In order to estimate the density profile according to Eq. (3), we fitted the observed visibility profiles with a power-law profile for uv distances ranging between 20 and 600 kλ for all coresexcept for core 1 in ISOSS J23053+5953. Here, we fitted the visibility profile between 60 and 600 kλ not considering the flat distribution toward the smaller baselines. Since the visibility profiles of cores 3 and 4 ISOSS J22478+6357 are flat, we do not fit their profiles. The temperature power-law index q, required for the calculation of p, is derived inSect. 4.5 and varies between 0.1 and 1 (Table 5).

The results for α and p are summarized in Table 5. The density power-law index ranges between 1 and 2. A detailed discussion and comparison with the density profiles of 22 cores in the more evolved CORE sample (Gieser et al. 2021) is given in Sect. 6.1.

Fig. 2 Multi-wavelength overview of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The top panels show a Spitzer IRAC composite RGB image at 4.5 (blue), 5.8 (green), and 8.0 μm (red). Emission with S∕N < 5 above the sky background is clipped and shown in black. The middle panels show Spitzer MIPS data at 24 μm. The bottompanels show a composite RGB image of Herschel PACS observations at 70 (blue), 100 (green), and 160 μm (red). Emission with S∕N < 5 above the sky background is clipped and shown in black. The Herschel clumps studied in Ragan et al. (2012) are indicated by green circles and labeled in green. In all panels, the NOEMA continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The area outside of the NOEMA field-of-view is gray-shaded. All mm cores identified in Beuther et al. (2021) are labeled in pink. The synthesized beam of the NOEMA continuum data is shown in the bottom left corner. A linear spatial scaleof 0.1 pc is indicated by a white scale bar.

3.3 Spectral energy distribution

The cold clumps toward both regions are studied by Ragan et al. (2012) using Herschel PACS observations at 70, 100, and 160 μm. Both regions were also observed with the Spitzer space telescope at 3.6, 4.5, 5.8, and 8.0 μm using the IRAC instrument and at 24 μm using the MIPS instrument. An overview of these archival data sets in comparison with the NOEMA 1.3 mm emission is shown in Fig. 2.

In the ISOSS J22478+6357 region, all mm cores are embedded within the R ≈ 0.1 pc sized Herschel clumps. A clump with strong 70 μm emission is detected toward the southwest with no associated 1.3 mm continuum emission (“clump 2”). This source also shows strong emission in the Spitzer IRAC and MIPS data. This suggests a more evolved source where the cold dust envelope was already disrupted by the protostar. The strongest MIR source is around the mm cores 13 and 15. These two mm cores, with a projected separation of ~4200 au, might in reality be the remaining dust envelope that is currently being disrupted by the protostar. Hennemann et al. (2008) find that this YSO is an intermediate-mass protostar. Cores 4, 5, and 6 also have bright counterparts in the Spitzer IRAC data.

The spatial extent of the two Herschel clumps in ISOSS J23053+5953 is larger, R ≈ 0.2−0.3 pc, compared tothe clumps in ISOSS J22478+6357. While cores 2, 3, 6, 7, 13, 14 are embedded in the northeastern clump peaking toward core 2, cores 1, 4, 5, 8, 9, 10, and 11 are associated with the southwestern clump peaking toward cores 1 and 5. MIR emission is detected around the following cores in ISOSS J23053+5953: cores 1, 2, 6, 7, 9, and 14. The southwestern clump is elongated toward core 12, but no significant FIR emission is detected there. There is also no corresponding MIR counterpart. Located at the edge of the primary beam with enhanced noise, core 12 might be an artifact.

The temperature and bolometric luminosity of the cold clumps are derived by Ragan et al. (2012) with Herschel PACS observations by fitting the SED. Following the nomenclature of Ragan et al. (2012), clumps 1− 5 toward ISOSS J22478+6357 and clumps 2 and 3 toward ISOSS J23053+5953 were covered by our observations (Fig. 2). The five clumps in ISOSS J22478+6357 have a temperature ranging between 18 and 21 K and bolometric luminosity between 7 and 10L_⊙. In the northeastern (clump 3) and southwestern clump (clump 2) in ISOSS J23053+5953 the dust temperature is 21 and 22 K and the bolometric luminosity is 441 L_⊙ and 869 L_⊙, respectively. However, the Spitzer data show that there is an additional warmer component.

We therefore refitted the SED of the cores with clear MIR counterparts in order to derive a more reliable estimate of the bolometric luminosity L. There are sources detected in the Spitzer data with no mm counterpart; however, in the following analysis we only focused on the sources with mm counterparts. It should be noted that while the Spitzer IRAC and NOEMA observations have a comparable angular resolution, the angular resolution of the Spitzer MIPS and Herschel PACS data is lower (Sect. 2.2). Smoothing the data to a common resolution would smear out all core features; instead, we performed a conservative cross-matching of the Spitzer MIPS sources and Herschel PACS clumps with the mm cores, as explained below.

In order to derive fluxes of the sources in the Spitzer IRAC and MIPS data, we performed aperture photometry using the photutils package. Background subtraction was performed by clipping sources with emission >3σ and then the median was computed and subtracted from the data. The DAOFIND algorithm (Stetson 1987) was used to identify sources with emission >5σ. A circular aperture with a radius of 2.″4 and 7.″6 was used to calculate the flux of the extracted sources in the IRAC and MIPS data, respectively. We did not subtract the background with a background annulus due to the crowded fields in the IRAC data. Instead we subtracted the background using the median value in the region as described above. The Spitzer sources are cross-matched to the closest mm core, if the projected distance is <1″ and <4″ for the IRAC and MIPS data, respectively. Aperture correction was applied to the derived fluxes, with a correction factor of 1.215, 1.233, 1.366, 1.568, and 1.61 for the fluxes at 3.6, 4.5, 5.8, 8.0, and 24 μm, respectively (Reach et al. 2005; Engelbracht et al. 2007).

In order to reliably fit the SED, we require that the mm cores are detected in three or four IRAC bands and therefore fitted the SED of evolved protostars. As a second further constraint, the core must either be detected at 24 μm or be associated with a Herschel clump peaking toward the position of the mm core (Table 4). This ensures that the mm fluxcan be accurately fitted with a cold component. For example, in the case of cores 6 and 7 in ISOSS J23053+5953, the mm cores have strong Spitzer IRAC counterparts; however, due to the poor angular resolution in the FIR, the Herschel clump (“clump 3” in Fig. 2) is not resolved toward these cores and peaks toward core 2 instead. While the warm component can be estimated with the Spitzer IRAC fluxes, a second cold component based on the mm flux cannot be constrained with only one data point. On the other hand, around core 13 in ISOSS J22478+6357, the Herschel clump (“clump 4” in Fig. 2) is not resolved also showing significant emission toward core 1. In this case, the Herschel PACS data points are not included in the fitting, but the strong 24 μm MIPS flux provided an additional data point peaking toward core 13. In summary, we can construct the SED from MIR to mm wavelengths for cores 4, 5, 6, and 13 in ISOSS J22478+6357 and for cores 1, 2, and 9 in ISOSS J23053+5953.

The fluxes of the Herschel PACS clumps are taken from Ragan et al. (2012), for which we cross-match the studied clumps to the closest mm core if one is detected (Table 4). The FIR observations resolve linear scales of 0.1−0.2 pc (Fig. 2), so the fluxes are integrated over a much larger area compared to the NOEMA 1.3 mm and Spitzer IRAC data. The 1.3 mm fluxes are taken from Beuther et al. (2021) and are listed in Table 4. The MIR and FIR photometricdata points are listed in Table 6 and the core SEDs are shown in Fig. 4.

For each core, we fit the SED with two or three components of a modified black body, $$B_{\nu }(T)=\frac{F_{\nu }\eta d^2}{{\kappa }_{\nu }M},$$(5)

in order to estimate the temperature and bolometric luminosity (see also Beuther & Steinacker 2007; Beuther et al. 2010; Linz et al. 2010; Ragan et al. 2012). The wavelength-dependent dust opacities, κ_ν, are taken from Ossenkopf & Henning (1994) for densities of 10⁶ cm⁻³ and 10⁵ yr of coagulation and with thin ice mantles. In consistency with all previous CORE and CORE-extension studies, we assume a gas-to-dust mass ratio of η = 150 ($\eta =\frac{M_{\text{gas}}}{M_{\text{H}}}$/$\frac{M_{\text{dust}}}{M_{\text{H}}}$, with $\frac{M_{\text{gas}}}{M_{\text{H}}}=1.4$ and $\frac{M_{\text{dust}}}{M_{\text{H}}}=0.0091$; Tables 1.4 and 23.1 in Draine 2011, respectively). This is higher than the typically assumed ratio of 100; however, recent observations suggest that η increases with galactocentric distance d_gal (Giannetti et al. 2017). These authors find a relation, where η = 130−320 at d_gal = 10 kpc. For most cores, two components were sufficient to model their SED; however, for core 2 in ISOSS J23053+5953, three components were required to reliably model all photometric data points.

The results of the temperature T_i and luminosity L_i of each i-th component, and total bolometric luminosity L are summarized in Table 6. The fitted total SED as well as each component are shown in Fig. 4. The contribution of cold dust with T ≈ 20−40 K is not sufficient to explain the observed fluxes at MIR wavelengths. In addition to the cold component, one or two warmer components with T = 400−600 K are required to reliably fit the SED (Table 6). This gives further evidence that some of the mm cores contain YSOs in a more evolved stage.

As an example for mm core 1 in ISOSS J23053+5953, two components with T = 23 K and T = 620 K are required to properly reproduce the FIR+mm and MIR fluxes, respectively. A strong 24 μm source is detected toward mm core 2 in ISOSS J23053+5953, for which three components with T = 28 K, T = 109 K, and T = 479 K are needed to describe the full SED.

One uncertainty arises from the fact that, even though the MIR and FIR sources are associated with the mm cores, the angular resolution of the observations from MIR to mm wavelengths vary and thus the fluxes are integrated over different angular sizes. This effect is particularly strong for the 160 μm fluxes whose observed fluxes are usually higher than the fluxes derived from the SED fit (Fig. 4).

Fig. 3 Radial visibility profiles of cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953. The azimuthally averaged visibility amplitudes $\overline{V}$, considering the real and imaginary components, are shown in black. The bottom x-axis shows the uv distance and the top x-axis is converted to the corresponding linear spatial scale. A fit to the data and its uncertainties (±1σ) are shown by the red solid and dashed lines, respectively.

4 Spectral line emission

In this section, we study the composition of the molecular gas by analyzing the 1.3 mm spectral line data. In the analysis of the original CORE data set (Gieser et al. 2021), the low-resolution data were used, which provided a continuous spectrum. However, as ISOSS J22478+6357 and ISOSS J23053+5953 are much colder and younger with typical line widths <3 km s⁻¹ (compared to a mean line width of 6 km s⁻¹ in the original CORE sample, Gieser et al. 2021), we used the high-resolution (~0.5 km s⁻¹) data if a high-resolution unit was placed toward the line. The properties of the spectral line data products (such as velocity resolution, synthesized beam, line noise σ_line) are summarized in Table 3.

Fig. 4 SED of mm cores with associated Spitzer IRAC emission in at least three bands. The fluxes (black points) are obtained from Spitzer IRAC (3.6, 4.5, 5.8, and 8.0 μm), Spitzer MIPS (24 μm), Herschel PACS (70, 100, and 160 μm), and NOEMA 1.3 mm observations (Tables 4 and 6). The total SED fit is shown by a red line and individual components are indicated by blue lines.

Fig. 5 CO 2−1 outflows in ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The blue and red contours show the integrated intensity of the blue- and red-shifted line wings, respectively. The intensity is integrated over vLSR − 21 to vLSR − 6 and over vLSR + 6 to vLSR + 21 km s−1 for the blue- and red-shifted integrated intensity, respectively. The corresponding velocity ranges are shown in each panel. Contour levels range from 10 to 90% of the peak integrated intensity with steps of 10%. In ISOSS J22478+6357 the red- and blue-shifted peak integrated intensity is 45 and 47 K km s−1 and in ISOSS J23053+5953 it is 112 and 140 K km s−1, respectively. Red- and blue-shifted outflow directions are indicated by red and blue arrows, respectively. The NOEMA 1.3 mm continuum data are shown in black contours. The dotted black contour marks the − 5σcont level. Solid black contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in black. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). Jet-like structures, seen in SiO 5− 4 emission, are indicated by orange ellipses (Fig. 7). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated bya black scale bar.

4.1 Molecular outflows

Molecular outflows are ubiquitous in low- and high-mass star-forming regions indicating indirectly the presence of gas accretion toward protostars (e.g., Beuther et al. 2002b; Wu et al. 2004; Zhang et al. 2005; Kölligan & Kuiper 2018). In low-mass star-formingregions, both the disk and outflow are commonly observed around protostars. Toward high-mass protostars, large bipolar outflows are found (e.g., Beuther et al. 2002b; Arce et al. 2007; Frank et al. 2014), but disk structures remain challenging to observe (e.g., Cesaroni et al. 2017; Ahmadi et al. 2019; Maud et al. 2019; Beltrán 2020; Johnston et al. 2020).

We used the integrated CO 2−1 emission in the blue- and red-shifted line wings to search for molecular outflows in ISOSS J22478+6357 and ISOSS J23053+5953. The blue- and red-shifted integrated intensity was computed from v_LSR − 21 km s⁻¹ to v_LSR − 6 km s⁻¹ and from v_LSR + 6 km s⁻¹ to v_LSR + 21 km s⁻¹, respectively.The results are shown in Fig. 5 for ISOSS J22478+6357 and ISOSS J23053+5953.

Both regions show multiple bipolar outflow signatures, indicated by arrows in Fig. 5. In Table 4 we list all mm cores with molecular outflows seen in CO 2− 1. We find outflow signatures toward cores 1−5 in ISOSS J22478+6357. The outflow toward core 1 is collimated. The outflow around core 2 has a quadrupolar morphology, similar to the outflow of the intermediate-mass YSO IRAS 22198+6336 (Sánchez-Monge et al. 2010). Either the quadrupolar shape arises from the cavity walls of a wide-angle outflow or two distinct outflows of unresolved multiple protostars are present. Toward the northwest of core 5 there could be a bipolar outflow with no associated mm continuum core. All outflow directions show a preferred orientation from the northeast to the southwest and are almost parallel with respectto each other. The orientation of the outflows are perpendicular to the filamentary gas traced by DCO⁺ for ISOSS J22478+6357 (Fig. A.2). Core 11 might have a north-south bipolar outflow, but the signature is not clear.

In ISOSS J23053+5953 a large-scale bipolar outflow is seen toward core 1 with the red- and blue-shifted lobes directed toward the northwest and southeast, respectively. A small projected outflow is observed toward core 9, which might be either a very young outflow just being launched or the result of an inclination effect. The CO 2− 1 line profile of core 9 is not significantly broader than the one of core 1, so an inclination effect is unlikely. Core 10 also hosts an outflow that is partially overlapping with the red-shifted outflow of core 1. Toward the northeast of ISOSS J23053+5953, red- and blue-shifted emission is detected around the location of cores 2, 3, 6, 7, 13, and 14; however, no clear outflow directions can be identified. For example, toward both cores 7 and 6 there is a weak signature of northeast-southwest bipolar outflows. We discuss in Sect. 4.3 that the region shows a steep velocity gradient toward these cores caused by a colliding flow, which is also seen in the CO line wing emission. It is therefore difficult to determine if these cores host outflows or not.

Clear detections of line wings from − 60 to − 40 km s⁻¹ were found toward ISOSS J23053+5953 using the CO 1−0 line observed with the IRAM 30 m telescope at an angular resolution of 21″ (Wouterloot et al. 1989). These authors derived the outflow properties and estimated a total outflow mass of M_out = 14.2−34.2M_⊙, momentum of P_out = 71.8−172.4M_⊙ km s⁻¹, energy of E_out = (40.5−97.4) × 10⁴⁴ ergs, a size of R_out = 1 pc for both blue- and red-shifted lobes, and an outflow timescale of τ_out = 10⁵ yr. Their outflow can be assigned to the large-scale outflow from mm core 1 with a red-shifted lobe directed toward the northwest and the blue-shifted lobe directed toward the southeast (Fig. 6 in Wouterloot et al. 1989).

Comparing the mass outflow rate of core 1, ${\dot{M}}_{\text{out}}=\frac{M_{\text{out}}}{{\tau }_{\text{out}}}\approx 2\times {10}^{-4}$ M_⊙ yr⁻¹ with the mass outflow rates of low- to high-mass YSOs (Fig. 11 in Henning et al. 2000), suggests that core 1 will form an intermediate- to high-mass star (see also Beuther et al. 2002b; Wu et al. 2004, 2005; López-Sepulcre et al. 2009; Maud et al. 2015b, for a comparison of the mass outflow rate and bolometric luminosity). The outflow has also been tentatively detected by Birkmann et al. (2007) using Plateau de Bure Interferometer (PdBI) observationsof CO 2−1 (we note that the red- and blue-shifted outflow directions are swapped in their Fig. 3).

The molecular outflows observed with NOEMA in both regions can be assigned to mm cores; however, their extent is larger than the observed mosaic, so a larger FOV would be required to properly derive the outflow properties (such as mass, momentum, energy, and timescale) of individual outflows. The multiplicity of bipolar outflows suggests that ongoing clustered star formation isoccurring in both regions. A recent study of the properties of the molecular outflow of core 1 in ISOSS J23053+5953 is presented in Rodríguez et al. (2021) using observations with the Submillimeter Array (SMA) and the VLA. Theseauthors detect compact centimeter (cm) emission toward the mm core, suggesting the presence of an ionized jet. The larger FOV of the SMA primary beam allowed them to properly derive the outflow properties, with M_out = 45.2 ± 12.6M_⊙ and τ_out = 1.5−7.2 × 10⁴ yr. The resulting outflow rate Ṁ_out ≈ 10⁻³ M_⊙ yr⁻¹ is even higher than previously determined by single-dish observations (Wouterloot et al. 1989). A further discussion of the protostellar outflows in ISOSS J22478+6357 and ISOSS J23053+5953 is given in Sect. 6.1.

Fig. 6 Line integrated intensity map of H2CO 30,3 − 20,2 of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The integrated intensity of the transition with a threshold of S∕N ≥ 3 is presented in color scale. The line intensity is integrated from a velocity of vLSR − 3 km s−1 to vLSR + 3 km s−1. It is marked if the transition is detected (√) or not (×), further explained in Sect. 4.2.1. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). Jet-like structures, seen in SiO 5− 4 emission, are indicated by orange ellipses (Fig. 7). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated bya black scale bar.

4.2 Spatial morphology and correlations

The spectral setup covers in total 26 emission lines detected in at least one of the two regions and an overview of the transition properties is shown in Table 2. We detect simple species consisting of 2− 5 atoms (CO, ¹³CO, C¹⁸O, ¹³CS, SO, H₂S, OCS, SO₂, H₂CO, H${}_2^{13}$CO, HNCO, HC₃N); deuterated ions (N₂ D⁺, DCO⁺); the shock tracer SiO; the cyclic molecule c-C₃H₂; and two complex organic molecules (CH₃OH, CH₃CN).

4.2.1 Line integrated intensity

In order to investigate the spatial morphology of each emission line, we computed the integrated intensity around the region systemic velocity from v_LSR − 3 km s⁻¹ to v_LSR + 3 km s⁻¹. The systemic velocity v_LSR for both regions is listed in Table 1. This velocity range covers 3 and 13 channels in the low- and high-resolution data, respectively. The integrated intensity (“moment 0”) maps for strong transitions of three key species, H₂CO (3_0,3−2_0,2), SiO (5− 4), and CH₃OH (4_2,3−3_1,2E) are shown in Figs. 6–8, respectively. The moment 0 maps of the remaining lines are shown in Appendix A. The noise in the integrated intensity maps σ_{int.intensity} was calculated by considering the velocity resolution (δv), line noise (σ_line), and number of channels (#channels): ${\sigma }_{\text{int}.\text{intensity}}=\delta v\times {\sigma }_{\text{line}}\times {(\#\text{channels})}^{0.5}$. In all integrated intensity maps, only locations with an integrated intensity > 3σ_{int.intensity} are presented.

The properties of the line data products, such as synthesized beam and line noise, are summarized in Table 3. It can be clearly seen in the integrated intensity maps that the noise is not uniform throughout the mosaic increasing toward the edge. For a comparison between the two regions, regardless of whether a transition shows significant emission throughout the FOV, we carefully investigated each line integrated intensity map, especially for the fainter lines, and searched for spatially resolved emission with S∕N > 3 that is not caused by noise artifacts. In Table 3 and in the integrated intensity maps we indicate for each region if the line is considered as detected (✓) or not (✗) in the integrated intensity map. For example, the OCS 18− 17 line (Fig. A.9) is clearly detected in ISOSS J23053+5953 toward core 1, 2 and 6, whereas in ISOSS J22478+6357 the line emission is irregular and caused by noise artifacts and therefore marked as not detected (✗). In Sect. 4.6, molecular column densities are derived from the core spectrum and there it was checked separately, if the transition has emission > 3σ_line.

While all transitions are detected in ISOSS J23053+5953, ISOSS J22478+6357 is more line-poor, with many transitions at higher upper energy levels (E_u/k_B ≳ 50 K) not detected (Table 3). In both regions, the emission of the three CO isotopologues (CO, ¹³CO, C¹⁸O, Figs. A.18, A.15, and A.10, respectively) is widespread across the FOV. The optically thick CO and ¹³CO lines trace the outer parts of the cloud and/or clump structure, while even the optically thin C¹⁸O emission is detected everywhere in the FOV. The detected H₂CO transitions (Figs. 6, A.7, and A.8) have extended emission in both regions. H₂CO is a good thermometer at temperatures <100 K and therefore its emission is used in Sect. 4.4 to create temperature maps of the regions.

In ISOSS J22478+6357, the line-richest object is core 1 with emission peaks of H₂CO, CH₃OH, SO, c-C₃H₂ (Figs. 6, 8, A.1, A.5, A.7, A.8, A.13, and A.17). The DCO⁺ and N₂ D⁺ emission (Figs. A.2 and A.21) peak toward core 2, but both molecules also have large-scale filamentary emission that is connecting the mm cores. H₂CO (Fig. 6) and SO (Figs. A.1 and A.13) also show emission peaks toward core 2. The cyclic molecule c-C₃H₂ has distinct emission peaks toward cores 1, 3, and 4 (Fig. A.5) tracing UV irradiated gas (e.g., Pety et al. 2005; Fontani et al. 2012; Mottram et al. 2020). There is no known PDR or UCHII region nearby and the emission peaks are co-spatial with the continuum peak of the cores. This suggests that the irradiation stems from the central protostar. Complementary observations at cm wavelengths would draw a clearer picture about the presence of any UCHII region.

We find locations in ISOSS J22478+6357 with strong molecular emission, but no 1.3 mm continuum counterpart. Toward the northeast and southwest of core 1 there are two molecular peaks (MPs), MP1 and MP2, seen clearly in H₂CO 3_0,3 −2_0,2 (Fig. 6) and SO 6₅−5₄ emission (Fig. A.13). The fact that these MPs are located at both sides of core 1 suggests that they are most likely shocked regions caused by the bipolar outflow of core 1 (Fig. 5). The dust emission toward the MPs can be too faint to be detected at our sensitivity limit, but due to shocks dust grain destruction of the mm -sized grains into undetectable μm-sized fragments might play a role as well. Toward the south of core 2, we also find a molecular emission peak only associated with faint 1.3 mm continuum emission (MP3). This is clearly seen in H₂CO 3_0,3 −2_0,2 (Fig. 6), CH₃OH 4_2,3 −3_1,2E (Fig. 8), and both SO transitions (Figs. A.1 and A.13). MP3 is connected to the red-shifted outflow cavity of core 2 (Fig. 5). ISOSS J22478+6357 is generally line-poor with simple species tracing the envelope, in which the cores are embedded. Molecular emission peaks with no mm continuum counterpart can be linked to molecular outflows, thus tracing shocked gas. Here, molecules, such as SO, H₂CO, and CH₃OH, which were initially frozen on the dust grains, are evaporated into the gas phase by the shock. Toward most of the cores, no distinct molecular emission is detected. This suggests that the cores are still too cold to have high molecular abundances in the gas phase.

ISOSS J23053+5953 is richer in line emission compared to ISOSS J22478+6357 (Table 3). Most of the molecular emission peaks at core 1, with some species peaking at cores 2 and 6. Core 1 has spatial emission peaks of SO (Figs. A.1 and A.13), ¹³CS (Fig. A.20), H₂S (Fig. A.3), OCS (Figs. A.9 and A.19), H₂CO (Figs. 6, A.7, and A.8), HC₃N (Figs. A.6 and A.23), and CH₃OH (Figs. 8, A.4, A.14, and A.17).

Core 2 has a prominent emission peak in the CH₃OH 5_1,4 − 4_2,3 line (Fig. A.4), but also in OCS (Figs. A.9 and A.19), HNCO (Fig. A.11), and CH₃CN (A.16) emission. Core 6 also shows many emission peaks and in comparison to core 1, strong emission of HNCO (Fig. A.11), c-C₃H₂ (Fig. A.5), and CH₃CN (Fig. A.16) is detected here, but also SO (Figs. A.1 and A.13), OCS (Figs. A.9 and A.19), and H₂S (Fig. A.3).

Similar to ISOSS J22478+6357, DCO⁺ is distributed throughout the FOV (Fig. A.2). Toward core 3, DCO⁺ emission peaks toward the north, while N₂ D⁺ emission (Fig. A.21) peaks toward the west. A possible scenario could be that there is a severe depletion of CO in the western position resulting in a low DCO⁺ abundance. This is further discussed in Sect. 6.2.

The cavity of the red-shifted outflow lobe of core 1 (Sect. 4.1) is seen between cores 1 and 5 in H₂CO (Fig. 6), CH₃OH (Figs. 8, A.4, A.14, and A.17), SiO (Fig. 7), and SO (Figs. A.1 and A.13), and HNCO (Fig. A.11) emission. The bipolar outflow is also seen as a dark lane in CO (Fig. A.18). In SiO (Fig. 7) and SO (Figs. A.1 and A.13), three jets can be identified toward the north of the region (labeled as Jet1, Jet2, and Jet3 in the integrated intensity maps). Jet1 is also seen in H₂CO emission (Fig. 6). While we do not find clear outflow signatures of the cores in this location, these jet features might indeed be caused by protostellar outflows (Sect. 4.1).

Toward the south of core 6 and toward the north of core 2, two emission peaks can be identified (MP4 and MP5) seen in H₂CO (Figs. 6, A.7, and A.8), H${}_2^{13}$CO (Fig. A.12), CH₃OH (Figs. 8, A.4, A.14, and A.17), SiO (Fig.7), SO (Figs. A.1 and A.13), SO₂ (Fig. A.22), HC₃N (Figs. A.6 and A.23), OCS (Figs. A.9 and A.19), and HNCO (Fig. A.11) emission. DCO⁺ emission is only enhanced toward MP5 (Fig. A.2). It could be that MP4 and MP5 are connected to potential molecular outflows as it is the case for MP1, MP2, and MP3 seen in SO, SiO, and H₂CO emission toward ISOSS J22478+6357. Unfortunately, the CO 2−1 line wing emission toward the location around cores 2, 3, and 6 is very complex and it is not possible to identify clear bipolar outflow signatures (Fig. 5). While we cannot rule out the presence of protostellar outflows causing these shocked regions, it coincides with a steep large-scale velocity gradient, for example, seen in NH₃ (Bihr et al. 2015) and DCO⁺ (Beuther et al. 2021), which hints at the presence of a colliding flow. This velocity gradient is further investigated in Sect. 4.3.

The emission of MP4 is elongated to the northeast toward mm core 6. All species having an enhanced emission toward MP4 and MP5 are known to trace shocked regions caused by protostellar outflows (Leurini et al. 2011; Benedettini et al. 2013; Moscadelli et al. 2013; Shimajiri et al. 2015; Palau et al. 2017; Tychoniec et al. 2019; Okoda et al. 2020; Taquet et al. 2020). There is enhanced c-C₃H₂ emission toward the north and south of core 2 (Fig. A.5) that could potentially trace, in addition to the colliding flow, a bipolar outflow.

A ring-like structure is seen in H₂CO emission in (Fig. 6) and in CO absorption (Fig. A.18) connecting the core 6, MP4, MP5, and core 2. This dynamically complex region might be an overlap of a colliding flow and the presence of protostellar outflows.

Fig. 7 Same as Fig. 6, but for SiO 5−4.

Fig. 8 Same as Fig. 6, but for CH3OH 42,3 − 31,2E.

4.2.2 Spatial correlations

To quantify the spatial correlation of the detected molecular emission lines, we applied the histogram of oriented gradients (HOG) method, for which a detailed description is given in Soler et al. (2019). In summary, HOG computes the relative orientation of the local intensity gradients of two position-position-velocity (PPV) cubes $I_{ijl}^{\text{A}}$ and $I_{ijm}^{\text{B}}$, where i and j run over the spatial axes and l and m over the spectral axes.

With the increasing bandwidth of correlators, sensitivity, and number of observed regions, it has become challenging to compare thespatial distribution of molecular emission and it is basically impossible to do that by eye on a channel-by-channel basis. The correlation function can also be used to study spatial correlations of the integrated intensity (e.g., Guzmán et al. 2018; Law et al. 2021). The HOG method also allows us to find a similar spatial morphology between two transitions. As the comparison is carried out in each velocity channel, potential velocity offsets between two molecules can be identified for example (Soler et al. 2019).

For the two regions in this study, we are able to compare the results obtained with HOG with the morphology of the integrated intensity maps (Sect. 4.2.1). However, for future line surveys toward star-forming regions and the analysis of the spatial morphology of the original CORE sample, HOG provides a convenient method to compare the line emission at high angular resolution and to find and study spatial correlations of molecular emission.

The projected Rayleigh statistic, (V), is a statistical test to determine whether the distribution of angles ϕ_ijlm between thegradients is nonuniform and peaked at a particular angle (see, e.g., Durand & Greenwood 1958; Batschelet 1972; Jow et al. 2018). In this application, the angle of interest is 0 degrees, which corresponds to the alignment of the iso-intensity contours in the PPV cubes.

In our application, we accounted for the statistical correlation brought in by the beam by introducing the statistical weights w_ij. Relative orientation angles of local intensity gradients ϕ_ijlm are computed after applying a Gaussian filter with kernel size Δ. The weighting is either $w_{ijlm}=\frac{\delta }{\Delta }$ with pixel size δ or w_ijlm = 0 in noisy regions. The projected Rayleigh statistic V_lm quantifies the amount of spatial correlation between the velocity channels l and m of two PPV cubes and is calculated from $$V_{lm}=\frac{{\sum }_{ij}{w}_{ijlm}\text{cos}(2{\varphi }_{ijlm})}{\sqrt{{\sum }_{ij}\frac{w_{ijlm}^2}{2}}}.$$(6)

Hence high values of V_lm correspond to a high spatial correlation, while low values of V_lm correspond to low or no spatial correlation, for example, when comparing two velocity channels dominated by noise. To recover the smallest spatial scales, we adopted a kernel size Δ equal to the synthesized beam. For both regions, we compare all pairs of detected emission lines listed in Table 3, except for the optically thick CO and ¹³CO transitions. For all transitions, the emission was compared in channels within a velocity range between v_LSR − 7.5 km s⁻¹ and v_LSR + 7.5 km s⁻¹. The only exception is the H₂CO 3_0,3 −2_0,2 line observed toward ISOSS J22478+6357, which lies at the edge of the high-resolution unit. Here we only considered a velocity range between v_LSR − 5.5 km s⁻¹ and v_LSR + 5.5 km s⁻¹. The peak projected Rayleigh statistic V for all transitions is shown in Fig. 9 for ISOSS J22478+6357 and ISOSS J23053+5953. A detailed example of two transitions with a high spatial correlation for each region is shown in Fig. 10.

In ISOSS J22478+6357, many transitions are not detected and thus excluded from the analysis (detections and non-detections for both regions are listed in Table 3). The low-energy transitions of H₂CO and SO have a high correlation. C¹⁸O, DCO⁺, and H₂CO trace large-scale emission. As expected, high correlations of transitions of the same species are found (H₂CO and SO). The high correlation between DCO⁺ and N₂ D⁺ is originating from the strong emission around mm core 2.

For ISOSS J23053+5953 we also find high correlations among combinations of transitions of the same species, which is the case for SO, H₂CO, and CH₃OH. While transitions with higher upper energy levels are generally less extended, the emission stems from the same location. The observed SiO transition shows a high correlation with both SO transitions toward the outflow of core 1, Jet1, and the shocked region caused by the putative colliding flow where MP4 and MP5 are located. With the exception of cores 1, 2, and 6 in ISOSS J23053+5953, the molecular emission does not peak toward the cores, which would be expected if the gas temperatures are already high enough to evaporate species frozen on the dust grains or allow efficient gas phase chemistry reactions. On the contrary, the emission of SO, SiO, H₂CO, and CH₃OH transitions peak toward the shocked region (MP4 and MP5).

In summary, with HOG we are able to find molecular species that show a high spatial correlation. This can be the case for chemically related species (e.g., H₂CO and CH₃OH) and for species tracing physical conditions such as a shock (e.g., SiO and SO). As expected, multiple transitions of the same molecule also have a high spatial correlation (e.g., for transitions of SO, H₂CO, and CH₃OH). In addition, the kinematic features of the molecular line emission can be studied in detail. We do not find significant velocity offsets for species with a high spatial correlation, but the two velocity components in ISOSS J23053+5953 can be clearly identified in Fig. 10 (right panel) at − 52 km s⁻¹ and at − 50 km s⁻¹ and are indicated by black dash-dotted circles in the figure.

Fig. 9 HOG results for ISOSS J22478+6357 (top panel) and ISOSS J23053+5953 (bottom panel). The peak projected Rayleigh statistic V is shown for each transition pair (Sect. 4.2.2). Examples of high spatial correlations in velocity space are shown in Fig. 10.

Fig. 10 HOG correlation results for DCO+ (3− 2) and H2CO (30,3 − 20,2) toward ISOSS J22478+6357 (top panel) and for SiO (5−4) and CH3OH (42,3 − 31,2E) toward ISOSS J23053+5953 (bottom panel). Emission maps of the channels with the highest correlation are shown in the left and middle panels. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated by a black scale bar. The projected Rayleigh statistic V, computed for all velocity pairs, is shown in the right panel. The red dashed line indicates equal velocity in both transitions. For ISOSS J23053+5953, the velocity regime around the two velocity components seen in line emission at v ≈ −52 km s−1 and − 50 km s−1 (Sect. 4.3) are highlighted by black dash-dotted circles.

4.3 Kinematic properties

A detailed study of the kinematic properties using the NOEMA + IRAM 30 m observations of ISOSS J22478+6357 and ISOSS J23053+5953 is presented in Beuther et al. (2021). DCO⁺ is a good tracer of the early stages of HMSF (e.g., Gerner et al. 2015).

Using the DCO⁺ 3−2 transition, both regions show distinct emission features (Fig. A.2). The integrated line intensity shows that the 15 cores in ISOSS J22478+6357 are connected by filamentary structures (Fig. A.2). Multiple velocity components are resolved toward substructures within the region and the line widths (FWHM) are small, on the order of ~1 km s⁻¹ (Beuther et al. 2021). The DCO⁺ line integrated intensity of ISOSS J23053+5953 does not show filamentary, but extended emission with many emission peaks close, but slightly offset from the core positions (Fig. A.2). There are two distinct velocity components, v ≈ −53 km s⁻¹ in the southeast direction and v ≈ −51 km s⁻¹ in the northwest direction (Beuther et al. 2021). This velocity gradient has already been reported by Bihr et al. (2015) using NH₃ emission at lower angular resolution and can be explained by a colliding gas flow triggering star formation. The line widths of individual DCO⁺ components are also small within the region, Δv ≤ 1.5 km s⁻¹, the only exceptions being cores 1 and 2 with Δv > 2 km s⁻¹ (Beuther et al. 2021). The thermal line width is 0.2 km s⁻¹ at 20 K (Beuther et al. 2021).

The presence of two velocity components is seen in various molecular tracers using single-dish (Wouterloot et al. 1988, 1993; Larionov et al. 1999)and interferometric observations (Bihr et al. 2015; Beuther et al. 2021). In Beuther et al. (2021) the DCO⁺ (3−2) intensity-weighted peak velocity (“moment 1”) map is shown to highlight the velocity gradients. For a more complete picture, we show in Fig. A.24 moment 1 maps of all observed lines with extended emission (¹³CO, C¹⁸O, SO, SiO, DCO⁺, H₂CO, and CH₃OH) with the exception of the optically thick CO 2−1 transition. The velocity gradient is clearly seen in all transitions. This velocity gradient is suggested to be caused by a colliding flow (Bihr et al. 2015; Beuther et al. 2021).

Another tracer of the dynamical processes in the regions is SiO. SiO is produced in shocked regions, for example, due to outflows, disk winds or converging gas flows. In a shock, silicon is sputtered off the grains and subsequently forms SiO in the gas phase (Schilke et al. 1997). Using the SiO 5−4 transition, significant emission is only detected toward the south in ISOSS J22478+6357 with no nearby continuum source. In Sect. 4.1, we find that this emission peak is shocked gas that is directed along the blue-shifted side of the outflow of core 1. In ISOSS J23053+5953, the spatial extent of the line integrated intensity of the SiO 5−4 is much larger peaking between cores 2, 3, and 6. Jet-like features (Jet1 and Jet2), which might be protostellar outflows originating from core 3, are seen toward the north of the region. SiO emission is also seen between cores 1 and 5 and likely caused by the outflow of core 1 (Sect. 4.1).

Fig. 11 H2CO parameter maps of ISOSS J22478+6357 derived with XCLASS. The fit parameter maps (Trot, N, Δv, voff) for H2CO are presented in color scale. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated by a black scale bar.

4.4 Formaldehyde distribution

The spectral line data can be used to derive molecular properties such as the rotation temperature T_rot and column density N in each pixel and to create parameter maps within the full FOV. It is computationally expensive to apply this method for all pixels and all detected molecules (Table 3); therefore, we only applied this pixel-by-pixel analysis to formaldehyde (H₂CO), for which we detect three transitions. For the regions of the original CORE sample, temperature maps are also derived using the high-density tracer CH₃CN (Gieser et al. 2021), but since it is not detected in ISOSS J22478+6357 and the emission is not extended in ISOSS J23053+5953, only H₂CO can be used to probe the gas temperature.

We employ the eXtended CASA Line Analysis Software Suite (XCLASS, Möller et al. 2017) to derive the following parameter set for a molecule: source size θ_source, rotation temperatureT_rot (“excitation temperature”), column density N, line width Δv, and offset velocity v_off with respect to the local standard of rest v_LSR. A detailed description of the used XCLASS setup is given in Appendix B.

The H₂CO rotation temperature maps are already presented in Beuther et al. (2021). Assuming T_kin ≈ T_rot, the H₂CO rotation temperature was used as input to estimate the H₂ column density N(H₂) and core mass M_core from the 1.3 mm continuum data (Table 4). Pixels with T_rot < 3 K are only found toward the edge of the FOV where the noise is high and the fits are unreliable.

The H₂CO parameter maps of ISOSS J22478+6357 are shown in Fig. 11. The rotation temperature T_rot is generally low varying between 10 and 50 K. Toward MP2 there is a ≈3″ × 3″ region with a high rotation temperature, T_rot > 150 K. As discussedin Sects. 4.1 and 4.2.1, this location is directed along the blue-shifted outflow of core 1 and is a shocked region with an enhanced H₂CO abundance and temperature increase. This can be seen in the column density map where N(H₂CO) is highest toward the core positions with strong mm continuum emission, but also toward this shocked region. The line widths are small, on the order of 1 km s⁻¹. The line width significantly increases at the positions of MP2 and MP3 with Δv > 3 km s⁻¹. A small east-west velocity gradient within the region is observed in the H₂CO envelope with v_off varying between ± 1 km s⁻¹.

The H₂CO parameter maps of ISOSS J23053+5953 are shown in Fig. 12. The average rotation temperature is slightly higher with an average of T_rot ≈ 70 K compared to ISOSS J22478+6357. A high increase with rotation temperatures > 150 K is observed toward MP4 and MP5, in cores 2 and 13, and toward cores 1 and 5. We note that the highest temperatures are not located toward the continuum peak positions of the cores, but are offset and/or between the cores. This implies that the cores themselves are still cold, in agreement with temperature estimates using NH₃ (Bihr et al. 2015), but shocks increase the local gas temperature and enhance the H₂CO abundance by a factor of two, and increase the line width. The shocked region in the northeast is associated with a colliding flow and potential outflows, while the shocked region between cores 1 and 5 is most likely associated with the strong outflow of core 1, which is also seen by strong velocity gradients in v_off (see Sect. 4.1 for a discussion of the outflows in the regions). Surrounding core 4 a ring-like structure of increase in T_rot, Δv, and v_off can be observed, which could be caused by the red-shifted outflow of core 1.

In summary, the H₂CO parameter maps reveal the large-scale gas distribution in both regions. We find that the temperature is enhanced toward the MPs that are caused by shocks. At this high angular resolution, we find that the temperature is higher in some locations compared to estimates using the Herschel PACS observations (~20 K, Ragan et al. 2012). Toward some mm cores and the shocked locations, the H₂CO column density and line width are enhanced as well. A small east-west velocity gradient is observed in ISOSS J22478+6357, while in ISOSS J22478+6357 there is a steep velocity gradient from the northeast to the southwest.

Fig. 12 Same as Fig. 11, but for ISOSS J23053+5953. Jet-like structures, seen in SiO 5−4 emission, are indicated by orange ellipses (Fig. 7).

Fig. 13 Radial temperature profile for cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953. The radial H2CO temperature profile is shown by the black data points. The outer radius rout, estimatedwith the clumpfind algorithm in the continuum data (Table 4), is indicated by the dash-dotted purple vertical line. A fit to the data and its uncertainties (±1σ) are shown by the red solid and dashed lines, respectively.

4.5 Radial temperature profiles

The derived H₂CO temperaturemaps can be used to create radial temperature profiles of the mm cores. To be consistent with the analysis of the density profiles (Sect. 3.2), we only applied the following analysis to cores with S∕N > 20. Radial temperature profiles were extracted from the H₂CO temperaturemap for cores 1−4 in ISOSS J22478+6357 and for cores 1−6 in ISOSS J23053+5953. We binned the data with a bin size 0.5× Δr corresponding to half of the average beam size ($\Delta r=\frac{{\theta }_{\text{maj}}+{\theta }_{\text{min}}}{2}$). The radial temperatureprofile for cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953 are shown in Fig. 13. Cores 3 and 4 in ISOSS J22478+6357 are close so their temperature profiles overlap.

We fitted the data assuming a power-law profile according to Eq. (4). We excluded the inner-most data point, which is diluted by the beam, and only fitted the data up to the bin that contains the outer radius r_out (Table 4). The outer radius r_out of each coreis taken from Beuther et al. (2021) and estimated by applying the clumpfind algorithm (Williams et al. 1994) to the 1.3 mm continuum data. For the characteristic radius and temperature we used r₅₀₀ = 500 au and T₅₀₀ = T(r₅₀₀). Flat temperature profiles (core 4 in ISOSS J22478+6357 and core 4 in ISOSS J23053+5953) or profiles that have an increasing temperature at increasing radii (core 2 in ISOSS J22478+6357 and core 3 in ISOSS J23053+5953) were not fitted.

The results for T(r₅₀₀), and q are summarized in Table 5. The temperature T(r₅₀₀) varies between 60 and 250 K indicating that some protostars have already heated up their surrounding envelope. The temperature power-law index q varies between 0.1 and 1.0. The results are consistent with results obtained in low- and high-mass cores (e.g., Shirley et al. 2000; Palau et al. 2014; Gieser et al. 2021). A comparison of the temperature profile with the cores of the CORE sample is discussed in Sect. 6.1.

For cores 1 and 2 in ISOSS J22478+6357, there is a significant temperature increase for r > r_out. The H₂CO temperature map (Fig. 11) reveals that toward the shocked gas of MP2 and MP3 the obtained rotation temperature is higher (>100 K) compared to the ambient gas (≈20−30 K). Thus, the radial temperature increase for r > r_out is related to the bipolar outflows of cores 1 and 2 (Sect. 4.1). The radial temperature increase of core 2 is also connected to the fact that it is located at the edge of the temperature map, where the temperature fits are unreliably due to an increase in noise (Fig. 11).

Flat temperature profiles are observed if the temperature gradient is not resolved. In addition, a flat profile can be observed if protostellar heating is not yet high enough to heat up the envelope. A high optical depth of the H₂CO lines could have an impact toward the densest regions of the cores, but since we can clearly derive temperature gradients for the densest cores (e.g., core 1 in ISOSS J23053+5953), line optical depth effects causing unreliable temperature estimates should not be an issue for the less dense cores. In Gieser et al. (2021) it was found that only in evolved HMSFRs containing HMCs the H₂CO temperature maps are not reliable in the densest central regions. In these cases, CH₃CN was used as a thermometer instead.

4.6 Molecular column densities

We quantify the molecular content by analyzing the molecular column densities N with XCLASS detected toward spectra of the 29 cores. The general setup of our XCLASS method is described in Appendix B. Here, we use the myXCLASSFit function to estimate the molecular column densities of detected molecules (Table 3). For each core, all species were fitted with one emission component. The column densities of the species are an input for the physical-chemical modeling of the cores (Sect. 5). We excluded CO and ¹³CO from the XCLASS fitting because of the high optical depth, and only fitted the optically thin C¹⁸O line. H${}_2^{13}$CO and SO₂ were also excluded, as there is no strong emission toward the 29 cores. Our physical-chemical model does not include shock chemistry; therefore, the column density of SiO was not derived either. It should be noted that the abundances of other species are also affected by shock chemistry (e.g., Palau et al. 2017), in the case of our observations we see shocked regions in enhanced SO, H₂CO, and CH₃OH emission (see Sect. 4.2.1).

Tables B.1 and B.2 summarize the best-fit column density for all 29 cores and species fitted with XCLASS. Upper limits are given, if the model peak intensity is <3σ_line, considering all transitions of a molecule. Upper limits are also given in cases where the upper column density error is larger than a factor of tenand when the lower column density error is zero. With these two additional constraints, unreliable fits can be further discarded.

As an example, we show the observed and XCLASS modeled spectra of core 1 in ISOSS J23053+5953 for all fitted transitions in Fig. B.2. Most of the transitions show a Gaussian-shaped line profile. The H₂CO 3_0,3−2_0,2 line is affected by self-absorption toward the systemic velocity and has broad line wings due to the bipolar outflow. The blue-shifted peak has a lower intensity than the red-shifted peak. This asymmetry arises from the outflowing gas being partially absorbed along the line-of-sight. The outflow is also seen in OCS 18−17 and 19−18 transitions having a weak double peaked line profile. The intensity of the CH₃OH lines that were only observed at a spectral resolution of 3 km s⁻¹ is overproduced by a factor of two in the XCLASS best-fit spectrum. The upper energy levels of these transitions are not significantly different from the remaining transitions observed with the high-resolution units. This mismatch is most likely due to the differences in spectral resolution, as the remaining two transitions are fitted well. In the low-resolution spectra, the emission of a few lines, for example for SO, is broader than the frequency range considered in the fit. This enabled that, in XCLASS, the non-Gaussian line wings caused by the outflows do not have a strong impact, which might cause broad features in the fitted line profile. For the remaining cores, the line widths are smaller and therefore the fit frequency range is sufficient to consider the emission of the full line. An example of the observed and XCLASS modeled spectra of core 1 in ISOSS J22478+6357 is presented in Fig. B.1, for which clearly fewer emission lines are detected, and hence column densities for less species can be estimated.

For all species, we computed the abundance relative to N(C¹⁸O) and abundance histograms are shown in Fig. 14. As the H₂ column density is derived from the NOEMA 1.3 mm continuum emission (Beuther et al. 2021) and suffers from spatial filtering, we use N(C¹⁸O) for a reliable comparison instead (see also Gieser et al. 2021). Each panel shows the results considering all cores, and separate histograms for ISOSS J22478+6357 and ISOSS J23053+5953 are presented as well. The logarithmic mean and standard variation of the abundance relative to N(C¹⁸O) considering both regions are shown in each panel. Some species are only detected toward ISOSS J23053+5953 (OCS, HC₃N, and CH₃CN). For SO, ¹³CS, DCO⁺, N₂ D⁺, H₂S, c-C₃H₂, and CH₃OH both regions have a similar abundance. However, the H₂CO abundance is about one order of magnitude higher in ISOSS J23053+5953 cores compared to the cores in ISOSS J22478+6357. A comparison of the mean abundance found toward the ISOSS regions compared to the CORE regions is given in Sect. 6.1.

5 Physical-chemical modeling

By analyzing the molecular line emission in Sect. 4.6, we find that the molecular content varies between the mm cores. Chemical timescales τ_chem of individual cores can be estimated with the physical-chemical model MUSCLE. The model was originally developed to study the evolutionary timescales of the IRDC, HMPO, HMC, and UCHII stages using single-dish observations of 59 HMSFRs(Gerner et al. 2014, 2015). The model consists of a static 1D physical structure coupled to the time-dependent gas-grain chemical network ALCHEMIC (Semenov et al. 2010). A detailed description of the physical-chemical model setupis given in Sect. 5.1. By definition, the chemical timescale of the initial IRDC stage is τ_chem = 0 yr, when the density reaches n = 10⁴ cm⁻³ (Gerner et al. 2014).

MUSCLE was applied to high angular resolution observations of the CORE pilot regions, NGC 7538S and NGC 7538IRS1 (Feng et al. 2016). In both regions, the cores have short chemical timescales of 11 000–12 000 yr. The model parameters were investigated in detail with CORE observations of the well-studied HMC AFGL 2591 VLA3 (Gieser et al. 2019), for which the derived chemical timescale is in agreement with estimates derived in the literature. A complete physical-chemical analysis of 22 cores in the CORE sample is presented in Gieser et al. (2021). The estimated chemical timescales varies from 20 000−100 000 yr. The chemical composition of the molecular gas of the cores varies even within a single region. An explanation for this variety is subsequent star formation. A few CORE targets contain UCHII regions, such as S106 or G139.9091, which both already show emission at cm wavelengths. Toward these regions, compact mm cores are found in the 1.3 mm continuum emission, but not many emission lines are detected, with the exception of CO isotopologues, SO, and H₂CO. In these cases, molecules are destroyed by the protostellar radiation and thus the spectra are line-poor compared to spectra toward HMCs of the CORE sample.

Fig. 14 Abundance histograms (relative to C18O) of all species fitted with XCLASS. The gray histogram shows the combined results of all 29 cores. The mean and standard deviation are shown in the top left corner. The green dashed and blue dash-dotted histograms show the results for the cores in ISOSS J22478+6357 and ISOSS J23053+5953, respectively.

5.1 Model setup

MUSCLE consists of a static spherically symmetric model core described by 40 radial grid points coupled the time-dependent gas-grain chemical network ALCHEMIC (Semenov et al. 2010). The density n(r) and temperature T(r) profiles are described by power-law profiles up to a radius r_out, according to Eqs. (2) and (4), respectively. From an inner radius r_in and inward, both density and temperature reach a constant plateau, n(r ≤ r_in) = n_in and T(r ≤ r_in) = T_in.

We used the same model parameters for the dust properties, interstellar radiation field, and initial chemical conditions used in the study of the 22 cores of the CORE sample, summarized in Table 5 in Gieser et al. (2021). Instead of running the model with only 100 logarithmic steps up to τ_model = 100 000 yr (Gieser et al. 2021), we used a finer grid of 999 logarithmic steps up to τ_model = 10⁷ yr. For the inner radius r_in we chose 500 au (Table 5). In Gieser et al. (2021), the outer radius r_out was determined from the temperature profiles, but for the ISOSS, we used the outer radius derived from the clumpfind analysis presented in Beuther et al. (2021) and summarized in Table 4.

We used the molecular hydrogen column density N(H₂) derived from the 1.3 mm continuum emission (Table 4). Short-spacing information is available for the spectral line data; however, the continuum data are obtained with NOEMA only, so missing flux can be an issue. Column densities derived with XCLASS of the following molecules were included in the model: C¹⁸O, ¹³CS, SO, OCS, DCO⁺, N₂ D⁺, H₂S, H₂CO, HNCO, HC₃N, c-C₃H₂, CH₃OH, and CH₃CN. The chemical network ALCHEMIC does not consider isotopologues except for a sophisticated deuterium (D) network, thus the column densities of C¹⁸O and ¹³CS were converted to the column density of their main isotopologue (CO and CS) using the isotopic ratio calculations from Wilson & Rood (1994) depending on the galactocentric distance d_gal. The ¹²C/¹³C isotopic ratio is 80 and 86 and the ¹⁶O/¹⁸O isotopic ratio is 607 and 649 for ISOSS J22478+6357 and ISOSS J23053+5953, respectively (Table 1). The ¹²C/¹³C isotopic ratios are in agreement with a study by Giannetti et al. (2014). These authors predict a ratio between 40−90 and between 45−95 at the corresponding galactocentric distance for ISOSS J22478+6357 and ISOSS J23053+5953, respectively.

A static spherically symmetric physical model is not sufficient to describe the protostellar evolution where the temperature and density are gradually increasing. However, coupling a chemical network with thousands of reactions to a hydrodynamical model is computationally expensive. As complex molecules, such as CH₃OH, are observed toward both regions, a simple chemical network would not be sufficient to describe the formation of all observed molecular species. Instead, we approximate the warm-up and density increase as step functions of several evolutionary stages with different temperature and density profiles. Using a sample of 59 HMFRS, Gerner et al. (2014, 2015) created template initial chemical conditions for the IRDC, HMPO, HMC and UCHII stage. For each core, we ran MUSCLE with IRDC, HMPO and HMC initial chemical conditions, referred to as the HMPO, HMC, and UCHII model, respectively, to derive the best-fit chemical timescale τ_chem. While we follow this nomenclature, it is not our goal to classify the cores into these evolutionary stages (see also Gieser et al. 2021).

For each initial condition model, the chemical age τ_chem was estimated from the best match of the modeled and observed column densities. The physical properties of each core were fixed and set to the observed physical properties (p, q, T₅₀₀, listed in Table 5). Therefore, the observed molecular column densities (Tables B.1 and B.2) constrained the best-fit model. As both regions are in general line-poor, only cores, for which at least eight of the in total 14 column density points are considered as detections and good fits and not upper limits, are modeled with MUSCLE. This includes four cores: core 1 in ISOSS J22478+6357 and cores 1, 2, and 6 in ISOSS J23053+5953.

Not all three initial condition models were considered for each core. The cores gradually warm up due to the protostellar radiation and a temperature decrease with increasing time would be unphysical. We therefore required that the observed temperature of the core at r_in is equal or higher than the temperature at r_in of the initial condition model, T_in ≥ T_in,init.

In each time step, the modeled radial abundance profiles are converted to beam convolved column densities with the same beam size as in the observations assumed. The modeled and observed column densities of the input molecules are compared by a least χ² analysis. We used a weighted χ², taking into account the percentage of well modeled molecules ϒ: ${\overline{\chi }}^2=\frac{{\chi }^2}{{{\rm Y}}^2}$. The model with the lowest ${\overline{\chi }}^2$ was then considered as the best-fit model. The column density of a molecule is well modeled if the modeled and observed column densities agree within a factor of ten. The chemical timescale is the sum of the timescale of the inital condition model τ_init and of the core model τ_model: τ_chem = τ_init + τ_model.

5.2 Chemical timescales

The best-fit chemical timescales are summarized in Table 5 for each modeled core. Detailed model results for each initial condition model are presented in Table C.1. Following the analysis of the mm cores of the original CORE sample (Gieser et al. 2021), we assume that the chemical timescales are uncertain by a factor of two (see also Gerner et al. 2014, 2015). A comparison of the observed and modeled column density is shown in Fig. 15 for all four modeled cores. Since in MUSCLE a factor of ten uncertainty for the modeled column densities is assumed and none for the observed column densities when deciding if a molecule is well modeled or not, we only show error bars for the modeled column density data points.

In the ISOSS J22478+6357 region, only core 1 can be modeled with MUSCLE and the estimated chemical timescale is 31 000 yr. This is in agreement with the fact that ISOSS J22478+6357 is young and still in a cold stage. With most molecules still frozen on the dust grains or not formed yet, the spectra appear line-poor. The modeled SO and H₂S column densities areboth higher than a factor of ten compared to the observed column densities. As discussed in Gerner et al. (2014) and in Gieser et al. (2021), the sulfur chemistry, especially SO, is poorly understood and is also affected by the presence of shocks. The elemental sulfur abundance had to be increased in the Gerner et al. (2014, 2015) models in order to properly explain the observed column densities of S-bearing species, but even in their case, the SO abundances could not be explained with MUSCLE. The same issue was found when modeling the 22 cores of the original CORE sample (Gieser et al. 2021). It is still under debate where and in what form sulfur is locked, for example, on the grains as OCS or H₂S or in the gas phase as neutral S (e.g., Goicoechea et al. 2021, and references within).

For core 1 in ISOSS J23053+5953, we estimate a chemical timescale of 19 000 yr. Both HMPO and UCHII model have a similar χ² of 0.458 and 0.457, respectively. However, the number of well-modeled molecules is higher for the HMPO model (71.4%) than for the UCHII model (64.3%), so the weighted ${\overline{\chi }}^2$ is lower for the HMPO model. In the physical-chemical analysis of the CORE sample, cores with τ_chem > 80 000 yr show emission from vibrationally excited lines that are not detected in the two ISOSS regions. In addition, the observed spectra toward the CORE cores with τ_chem > 80 000 show typicalhot-core type spectra and with more emission lines detected compared to core 1 in ISOSS J23053+5953. This suggests that the HMPO model is indeed the more realistic model. The best-fit chemical age for cores 2 and 6 are similar to the chemical timescale of core 1, 19 000 and 18 000 yr, respectively.

The MUSCLE model overestimates the SO column density for all cores except for core 6 in ISOSS J23053+5953. The H₂CO column density is underestimated by MUSCLE for all cores except core 1 in ISOSS J22478+6357. These species are difficult to model with MUSCLE since the sulfur chemistry was already an issue in the initial abundance models by Gerner et al. (2014, 2015) and shock chemistry is not included, which is important to consider, for example, for outflow sources having enhances SO and H₂CO emission (see also Sect. 4.2.1). A comparison of the chemical timescales τ_chem of the four CORE-extension cores with the 22 cores of the CORE sample is discussed in Sect. 6.1.

Fig. 15 Comparison of the observed and modeled column densities of all cores modeled with MUSCLE. The observed column density is shown in red and upper limits are indicated by arrows (Tables B.1 and B.2). The column densities of the MUSCLE best-fit model are shown by the blue data points assuming a factor of ten uncertainty.

6 Discussion

The core properties of both regions are discussed and compared with the results of the 22 cores of the original more evolved CORE sample (Gieser et al. 2021) in Sect. 6.1. A more detailed discussion of the complex dynamics in ISOSS J23053+5953 is given in Sect. 6.2.

6.1 Core properties

Both regions, located in the outer galaxy (d_gal ≈ 10 kpc), fragment into a total of 29 mm cores (with S∕N ≥ 5) as revealed by the high angular resolution observations with NOEMA at 1.3 mm (Beuther et al. 2021). Herschel observations at FIR wavelengths show that the mm cores are embedded in cold clumps (Ragan et al. 2012). MIR observations with the IRAC and MIPS instruments on the Spitzer space telescope reveal that multiple mm cores have MIR counterparts (Fig. 2). This suggests that in both regions star formation is occurring subsequently with some more evolved sources already becoming bright at MIR wavelengths. The strong mm cores with no or a weak associated MIR sources are likely younger. In this case, the luminosity of the protostars is still too faint and/or the protostars are still too embedded.

6.1.1 Clustered star formation

The total clump masses of the cold gas are 104 M_⊙ and 488 M_⊙ for ISOSS J22478+6357 and ISOSS J23053+5953 (Ragan et al. 2012). Tackenberg et al. (2012) estimate that a gas reservoir of 1000 M_⊙ is required to form a 20 M_⊙ star on clump scales at a size of R ≈ 0.3−1.0 pc. This suggests that in ISOSS J22478+6357 low- and intermediate-mass stars can form. This is in agreement with the results of Hennemann et al. (2008) who estimate that the bright MIR source toward mm core 13 is forming an intermediate-mass star. However, toward ISOSS J23053+5953, core 1 with its strong bipolar outflow (Wouterloot et al. 1989; Rodríguez et al. 2021) and being embedded within a large gas reservoir (Birkmann et al. 2007) is an accreting intermediate-mass protostar, destined to become a high-mass star (see, e.g., Beuther et al. 2007a).

The SED fitting (Sect. 3.3) shows that for individual cores the bolometric luminosity is not only dominated by the FIR emission of the cold dust, but a significant contribution can arise at MIR wavelengths. In ISOSS J22478+6357, the strong mm cores 1 and 2 have no associated MIR source, while in ISOSS J23053+5953 mm cores 1 and 2 can be associated with MIR emission. Some of the faint mm cores show MIR emission, for example, core 13 in ISOSS J22478+6357 is the brightest in the MIR. As both regions harbor young intermediate- to high-mass protostars (mm-bright) surrounded by more evolved low-mass YSOs (MIR-bright) suggests that clustered subsequent star formation occurs. It seems that first a generation of low-mass stars form, with bright emission in the near-infrared (NIR) and MIR, and then higher-mass protostars, bright at mm wavelengths, form afterward. As an example, in ISOSS J22478+6357 the MIR-bright source toward mm core 13 has an estimated age of ~ 1−6 Myr (Hennemann et al. 2008) and for the nearby mm core 1 we derive a shorter chemical timescale of 31 000 yr (Table 5). While the angular resolution of the Spitzer IRAC and NOEMA observations is comparable and resolve the core features, the lower angular resolution of the Spitzer MIPS and Herschel PACS data trace the clump scales, in which multiple cores can be embedded (Fig. 2). Unfortunately, high angular resolution observations at FIR wavelengths, which would significantly improve the SED fit and bolometric luminosity estimate of cold components toward individual mm cores, are currently not possible (e.g., Linz et al. 2021).

In both regions a small cluster is undergoing active star formation (Fig. 2). Kumar et al. (2003a) studied young star-forming regions that are classified as precursors to UCHII regions with no free-free emission detected at radio wavelengths. Using observations in the NIR and FIR, these authors found that the regions form stars in a clustered mode. Toward the young massive protocluster IRAS 22134+5834, a massive protostar shows strong emission in the FIR and has a weak radio continuum, while the surrounding protostars on the other hand are bright in the NIR (Kumar et al. 2003b; Palau et al. 2013; Wang et al. 2016). This is an indication for an early phase of massive star formation with protostars at a range of evolutionary stages. This is also observed toward both ISOSS regions in our study where strong mm cores and bright MIR sources coexist (Fig. 2). The formation and evolution of massive stars is much faster, but it is also commonly observed that high-mass protostars are surrounded by a cluster of more evolved low-mass protostars (see also e.g., Kumar et al. 2006).

This is in agreement with the scenario proposed by Kumar et al. (2020) that low-mass star formation occurs in filaments, while HMSF occurs in hubs, which take longer to assemble. Thus, high-mass protostars are observed to be surrounded by more evolved low-mass protostars. We find that in ISOSS J22478+6357, the mm cores are embedded within filamentary structures, whereas in ISOSS J23053+5953 the system consist of two bigger clumps resembling more of a hub (Sect. 4.2.1). This also supports the theory that star formation occurs at the intersection of filaments and gas hubs, and that the growing massive cores could be fed by fresh matter from and/or via these filaments (e.g., Myers 2009; Peretto et al. 2013; Tigé et al. 2017; Kumar et al. 2020). However, it is difficult from an observational point of view to observe this possible gas accretion from large to small scales. Both ISOSS J22478+6357 and ISOSS J23053+5953 seem to be embedded in a filament-hub system as observed with the Herschel Spectral and Photometric Imaging REceiver (SPIRE). But currently no velocity information of these structures is available and thus follow-up line observations of a larger FOV are required in order to investigate the clustered star-formation properties in both regions.

Another explanation could be that additional external feedback triggers a further generation of stars. In ISOSS J23053+5953, the colliding flow could have triggered further star formation at the colliding front toward cores 2, 3, and 6. In ISOSS J22478+6357 the nearby IRAS source, which could be an expanding HII region, but has not been studied, could have an impact on the cluster.

6.1.2 Physical properties

In both regions, clear signatures of protostellar outflows are observed in line wings of the CO 2−1 transition (Fig. 5). For comparison, Zhang et al. (2007) find multiple collimated outflows toward the high-mass protocluster AFGL 5142 with CO, SO, and SiO observations. In AFGL 5142, the bipolar outflows are randomly oriented and the outflows are associated with mm cores. In the filamentary IRDC G28.34+0.06 P1, Wang et al. (2011) detect mm cores aligned in a filamentary structure. Even though the cores show no line emission except for CO, all cores harbor bipolar jet-like CO outflows. The outflows are oriented almost parallel to each other and perpendicular to the filament structure. The molecular outflows are observed in SO, SiO, H₂CO, and CH₃OH emission as well (Zhang et al. 2015).

This is similar to the line-poor mm cores toward ISOSS J22478+6357, for which the outflows are aligned and perpendicular to the filamentary gas as well (Fig. 5). This suggests that for protostars forming along filaments, there is a preferred orientation of the outflows and that the filaments can provide gas that can be accreted. The magnetic field morphology is a topic for future investigations since the magnetic field is typically also oriented perpendicular to the high-density filaments (e.g., Soler et al. 2017).

On the other hand, the outflows that can clearly be identified in ISOSS J23053+5953 do not have a preferred orientation. There is no clear signature of a filament, but rather the presence of two massive clumps (Fig. 2) and a colliding flow toward the location of the mm cores (Fig. A.24). Similar to AFGL 5142, outflows toward more massive star-formingregions are influenced by the gravitational interaction between the protostars, which might cause the outflows to precess and change their orientation with time (for detailed simulations of high-mass outflows, see, e.g., Peters et al. 2014).

The physical structure, specifically the radial temperature and density profile, is analyzed for the ten strongest cores within both regions. The 1.3 mm visibility profiles (Fig. 3) were used to estimate the density power-law index p of five cores(Table 5). We derive a mean density power-law index of p = 1.6 ± 0.4. The same analysis was applied to 22 cores of the original CORE sample (Gieser et al. 2021), with a mean of p = 2.0 ± 0.2. This suggests that from the IRDC to HMC stage, the density profile does not change significantly within these short evolutionary timescales. In Gieser et al. (2021) the mean density power-law index was compared to results at various spatial scales in the literature. It is found that p does not change from clump to core scales. A mean value of p ≈ 2 is also found toward low-mass star-forming regions (Motte & André 2001). Hence a density profile around two appears to be present over a broad range in evolutionary timescale as well as protostellar mass. The density profile of an initial (high-mass) starless core and in a more evolved UCHII region phase might differ from the intermediate IRDC, HMPO, and HMC stage, though.

The radial temperature profile of the cores (Fig. 13) is derived from the H₂CO rotation temperature maps (Figs. 11 and 12). The mean temperature power-law index is q = 0.4 ± 0.3 ranging between q = 0.1−1 within the six mm cores with a radially decreasing temperature profile (Table 5). In some cores, we do not observe a radial temperature gradient, which might not be resolved by our observations. For the cores of the original CORE sample, H₂CO and, in addition, CH₃CN temperature maps were used to estimate the temperature profiles of individual cores. The mean temperature power-law index of the 22 cores of the original CORE sample is q = 0.4 ± 0.1, which is predicted from theoretical calculations (Emerson 1988; Osorio et al. 1999; van der Tak et al. 2000) and similar to the value found toward the ISOSS cores.

It is expected that the initially cold cores develop a radially decreasing temperature profile due to protostellar heating. The fact that the density and temperature power-law index are similar from the IRDC to HMC phase suggests that the density and temperature gradient form very early on and stay roughly constant. Since observational constraints of the temperature and density profiles at core scales are still scarce, our findings provide useful constraints for physical models and simulations, but in addition, the density and temperature structure are crucial for chemical models as well (Sect. 5).

Comparing results of the H₂CO temperature maps (Table 5) and the temperature of individual components in the SED fit (Table 6) reveals a temperature gradient of decreasing gas temperature at increasing spatial scale. The Herschel observations tracing the larger scale clumps at 0.1−0.3 pc reveal a cold component of ~20 K. Thus, the mm cores form within a cold clump. The CORE-extension observations show that at a spatial scale of ≈ 3000−4000 au the gas temperature increases toward the cores, with H₂CO rotation temperatures up to ~200 K. The MIR sources appear as point sources in the Spitzer IRAC observations. The strong MIR emission arises from the close vicinity of the protostar and reveals a hot component of 400−600 K toward the mm cores with MIR counterparts.

6.1.3 Molecular gas

Both ISOSS J22478+6357 and ISOSS J23053+5953 appear line-poor compared to many regions in the original CORE sample. In total, 26 transitions can be identified (Table 2), all of them toward ISOSS J23053+5953, while toward ISOSS J22478+6357 only 14 transitions are detected (Table 3). Toward the regions of the CORE sample, strong line emission of SO₂, vibrationally excited HC₃N and torsionally excited CH₃OH and less abundant isotopologues, such as O¹³CS and ³⁴SO₂ are commonly detected in the spectra. In addition, more transitions of HNCO and CH₃CN with higher upper energy levels are detected. This contrast is in agreement with the fact that the gas temperature and density in ISOSS J22478+6357 and ISOSS J23053+5953 are lower suggesting that these regions are less evolved.

The mm cores are embedded in an envelope seen in H₂CO and DCO⁺ emission (Figs. 6 and A.2). While both regions have similar clump properties (Ragan et al. 2012), the high angular resolution data reveal complex dynamics that affect the molecular gas of the regions (Sects. 4.2.1 and 4.3).

While ISOSS J22478+6357 is more quiescent and line-poor, a colliding flow in ISOSS J23053+5953 is creating a large shocked region that is seen in line-rich emission knots (MP4 and MP5), which are not associated with the mm cores and are in a steep velocity gradient (Fig. A.24). In both sources, the temperatures are low throughout the region, so most of the molecules around individual cores are either still frozen onto the dust grains or have not yet formed. Shocked regions are caused either by bipolar outflows (MP1, MP2, and MP3 in ISOSS J22478+6357, Fig. 5) or by a large-scale colliding flow (MP4 and MP5 in ISOSS J23053+5953, Fig. A.24). These shocks are traced by emission of SO, SiO, H₂CO, and CH₃OH (e.g., Figs. A.1, 7, 6, and A.4). The presence of molecular outflows, which mightalso create shocks, toward mm cores surrounding MP4 an MP5 cannot be ruled out.

The mean abundance, relative to N(C¹⁸O) and standard deviation of the 29 cores in both ISOSS regions in comparison with the mean abundance of the 22 cores in the CORE sample is shown in Fig. 16. Only abundances of species that were analyzed in both studies are shown. While there is clearly a scatter for the ISOSS and more evolved CORE cores, which is most likely caused by individual cores being at slightly different evolutionary stages, the mean abundances of all species of the ISOSS cores are about an order of magnitude lower than the mean abundances of the CORE cores. This can be interpreted in an evolutionary sense that most of the species are still depleted early on and the gas phase abundances and gas density increase during the evolution of the cores.

The molecular emission in the large mosaics allows us to investigate different physical processes. The larger scale molecular cloud is traced by CO isotopologues (CO, ¹³CO, and C¹⁸O, Figs. A.18, A.15, and A.10), with even C¹⁸O being detected everywhere in the FOV. We find that DCO⁺ and H₂CO are also good tracers of the larger scale emission, that is, the filamentary structure in ISOSS J22478+6357 and the clumps in ISOSS J23053+5953 (Figs. A.2 and 6). Distinct molecular emission peaking toward the mm cores (SO, ¹³CS, OCS, DCO⁺, N₂ D⁺, H₂S, H₂CO, H${}_2^{13}$CO, HNCO, HC₃N, c-C₃H₂, CH₃OH and CH₃CN) is seen for cores 1 and 2 in ISOSS J22478+6357 and for cores 1, 2, and 6 in ISOSS J23053+5953. Bipolar molecular outflows are traced by broad CO line wing emission (Fig. A.18), SO (Figs. A.1 and A.13), SiO (Fig. 7), H₂CO (Fig. 6), H${}_2^{13}$CO (Fig. A.12), and CH₃OH (Fig. 8).

Multiple molecular emission peaks in molecular emission are detected in both regions (MPs 1−5). Qiu & Zhang (2009) detect a MP in the HMSFR HH 80-81 also with no mm continuum counterpart. In their case, a CO outflow arises from this molecular core, while in our case, the molecular emission, which also is not associated with any mm continuum, stems from shocks caused by outflows and a colliding flow.

Fig. 16 Mean abundances (relative to C18O) and standard deviation of the cores in both ISOSS regions (green circles) in comparison with the cores in the CORE sample (blue triangles) taken from Gieser et al. (2021). Only species that were studied both in this work and in Gieser et al. (2021) are shown.

6.1.4 Physical-chemical modeling results

In Gieser et al. (2021), the physical-chemical model MUSCLE is applied to 22 cores in order to estimate the chemical timescale, τ_chem, and a mean of 60 000 yr, with a spread of 20 000–100 000 yr, is estimated. It is not possible to apply MUSCLE to all cores in ISOSS J22478+6357 and ISOSS J23053+5953 as the spectra are in general line-poor. The mean chemical timescale of the ISOSS cores is 20 000 yr. The 1.3 mm spectral setup is ideal to study the molecular emission around cores in the HMPO and HMC stage (Gieser et al. 2021), but in the IRDC stage spectra appear line-poor (Figs. B.1 and B.2). Complementary observations at 3 mm with many simple species, including their deuterated counterparts, in the ground state levels would help to better constrain the chemical timescales of the mm cores in ISOSS J22478+6357, but also of the line-poor mm cores in the original CORE sample.

For ISOSS J22478+6357 we are only able to model core 1 that shows the most emission lines compared to the remaining cores in the region. As expected from the low temperature in the region, the estimated chemical timescale is low (~31 000 yr). The remaining strong mm cores hosting bipolar outflows as well, are most likely in an early evolutionary stage, where more species are still frozen onto the grains, which is also seen in low kinetic temperatures of the cores (Table 4). The estimated chemical timescale of core 1 in ISOSS J22478+6357 is ≈ 10 000 yr higher compared to the mm cores in ISOSS J23053+5953. Since M_core and T₅₀₀ are low, this low-mass protostar cannot reach densities as high as in core 1 in ISOSS J23053+5953 in order to have a more complex molecular content. The ISOSS J22478+6357 cluster seems slightly older, with a MIR-bright YSO around mm core 13 at an estimated system age of 10⁶ yr (Hennemann et al. 2008) already present.

Core 1 in ISOSS J23053+5953 is most likely a young intermediate- to high-mass protostar with a strong bipolar outflow at an estimated chemical timescale of ≈20 000 yr. The single-dish observations reported in Vasyunina et al. (2014), centered on core 1, show the presence of a few more complex organic molecules, such as CH₃CCH and CH₃CHO, detected at 2 and 3 mm. This indicates the onset of some more complex chemistry toward this particular core. Unresolved CH₃CN emission in our NOEMA data toward the continuum peak already suggests that the inner parts already reach high temperatures and densities (Fig. A.16). For cores 2 and 6 in ISOSS J23053+5953, we also estimate a comparable chemical timescale. Both can be associated with MIR sources as well (Fig. 2), both significantly brighter in the MIR than core 1, and also have distinct MPs in the line integrated intensity maps (Sect. 4.2.1).

6.1.5 Correlations

Using the CORE and CORE-extension observations, the physical and chemical properties of 32 fragmented mm cores were studied in detail (ten cores in this work and 22 cores in Gieser et al. 2021). The following derived properties are investigated for correlations: the 1.3 mm peak intensity and integrated flux (I_peak and F₁₃₀₀), the kinetic temperature estimated from the H₂CO temperaturemap (T_kin), core mass M_core, molecular hydrogen column density N(H₂), kinetic temperature at r = 500 au obtained from the radial temperature profile T₅₀₀, visibility, density, and temperature power-law index (α, p, and q), chemical timescale (τ_chem), and M/L ratio.

The absolute flux calibration of the NOEMA data is expected to be 20% or better. We therefore assume a 20% uncertainty for I_peak and F₁₃₀₀. Using the MCMC error estimation algorithm, the uncertainty of the H₂CO rotation temperature estimated with XCLASS is ~30% (Gieser et al. 2021). Since T_kin is approximated by T_rot(H₂CO), we assume a 30% uncertainty for T_kin. The core mass and H₂ column density are uncertain by a factor of three (taking into account further uncertainties of, for example, the distance, gas-to-dust mass ratio, and dust absorption coefficient, Beuther et al. 2021). The uncertainties for T₅₀₀, α, p, and q are taken from the uncertainty of the fit (Table 5). The chemical timescales are estimated to be uncertain by a factor of two (Sect. 5.2). The M/L ratio (M/$L=\frac{M_{\text{core}}}{L}$, Table 6) is assumed to be uncertain by a factor of three. The M/L ratio of the CORE cores, based on single-dish observations and taken from Beuther et al. (2018), is an average of the whole region. But for the ISOSS cores we are able to constrain the M/L ratio on core scales (Sect. 3.3). For completeness, we also show the clump-averaged M/L ratio taken from Ragan et al. (2012) of the ISOSS cores that have an associated Herschel clump peaking toward the core (core-clump associations are listed in Table 4). The clump-averaged M/L ratio of the ISOSS cores is shown in red in Fig. 17, but not considered in the calculation of the correlation coefficient.

In order to find correlations or anticorrelations, we computed the Spearman correlation coefficient r_S (Cohen 1988) for all parameter combinations. If r_S = 1, a positive correlation exists and if r_S = −1 an anticorrelation exists between the parameters. For uncorrelated parameters, r_S = 0 holds. All parameter pairs with |r_S|≥ 0.4 are shown in Fig. 17. The ISOSS and CORE cores are shown in green and blue, respectively.

A high correlation exists between I_peak and F₁₃₀₀, but the figure reveal that the cores of the ISOSS regions are fainter than the cores of the CORE regions with a small overlap in between. This correlation spans over four orders of magnitude highlighting that the CORE-extension cores are a complement of the CORE cores. Brighter mm cores also have a higher kinetic temperature T_kin. As N(H₂) and M_core are computed from I_peak and F₁₃₀₀, respectively, we find a high correlation here as well. Thus, M_core and N(H₂) also have a high correlation. T_kin and T₅₀₀ correlate with N(H₂) and M_core, respectively.

The temperature index q is anticorrelated with the density index p, such that cores with flatter temperature profiles have a steeper density profile. A strong correlation between α and p can be attributed to Eq. (3). We do not find correlations for α, p, and q with other physical parameters. The mean indices for the density and temperature power-law indices are similar in both studies, implying that these parameters do not change significantly from the IRDC to HMC stage.

Both the peak intensity and H₂ column density of the mm cores correlate with τ_chem such that brighter mm cores with higher H₂ column densities have a higher estimated chemical timescale. The M/L ratio is a tracer of the evolutionary trend such that more evolved protostars have a lower M/L ratio (e.g., Sridharan et al. 2002; Molinari et al. 2008, 2010, 2016, 2019; Maud et al. 2015a; Urquhart et al. 2018). The peak intensity, integrated flux, and H₂ column density of the mm cores are anticorrelated with the M/L ratio. As would be expected as a tracer of the evolutionary trend, the chemical timescale τ_chem is also anticorrelated with the M/L ratio.

The correlation of the physical parameters of the mm cores in the CORE and CORE-extension sample reveal that some of the studied physical properties change with time: cores with a high mm brightness also have a higher temperature, core mass, H₂ column density and lower M/L ratio. This suggests that the initial low-mass cores gain mass with time while also increasing the gas temperature and density. Both the M/L ratio and chemical timescales τ_chem are tracers of the evolutionary sequence. It is difficult to determine the bolometric luminosity L of cores, which are clearly resolved in MIR and mm observations, but not at FIR wavelengths (Fig. 2). Since it is currently impossible to reach subarcsecond resolutions at FIR wavelengths (e.g., Linz et al. 2021), more reliable SED fits remain challenging, especially in crowded regions (Sect. 3.3). However, chemical timescales estimated with MUSCLE are determined by the observed molecular column densities, for which subarcsecond resolution observations are possible with current interferometers. We do not find that the indices of the temperature and density power-law profiles have a high correlation with the remaining parameters suggesting that they are set very early on and do not change significantly from the IRDC to HMC stage.

Fig. 17 Core parameters with a Spearman correlation coefficient of |rS|≥ 0.4. The results for the cores in ISOSS J22478+6357 and ISOSS J23053+5953 are shown in green (this work, Tables 4–6) and the results for the cores of the original CORE sample are presented in blue (Gieser et al. 2021). In each panel, the correlation coefficient rS and p-value are shown in the top left corner. For completeness, the clump-average M/L ratio taken from Ragan et al. (2012) of the corresponding ISOSS cores (Table 4) is shown in red.

6.2 Dynamics in ISOSS J23053+5953

In ISOSS J23053+5953 a colliding gas flow creates enhanced abundances of shock tracers such as SO, ¹³CS, SiO, OCS, SO₂, H₂CO, H${}_2^{13}$CO, HC₃N, CH₃OH, and CH₃CN (indicated as MP4 and MP5 in the integrated intensity maps, Sect. 4.2.1). These shocked locations reveal molecules that were initially depleted. Therefore, it allows us to investigate which molecules already resided on the grains. In addition, SiO is created by Si sputtered off the grains and forming SiO in the gas phase (Schilke et al. 1997). These emission peaks are spatially located at a velocity gradient seen clearly for all species with extended emission (Fig. A.24). However, since both cores 2 and 6 are young protostars, molecular outflows are expected to be present, but difficult to identify (Fig. 5).

Core 2, which is embedded within a common dust envelope with core 3 (Fig. 1), shows emission of c-C₃H₂ (Fig. A.5) toward the north and south of the mm core. This might trace an outflow cavity, but a clear outflow direction cannot be inferred due to the complex dynamics. We therefore conclude that a colliding flow creates a shocked region between cores 2, 3, and 6, but protostellar outflows are most likely present as well.

Colliding flows are found in other regions of the Galaxy. For example, in the well-known giant molecular cloud complex NGC 6334, colliding flows are reported in the source known as “V” by Juárez et al. (2017), where two velocity components separated about 2.4 km s⁻¹ are seen in H¹³CO⁺, with SiO emission at the intermediate velocities and also at the interface of the two flows (Fig. 7 of Juárez et al. 2017). Multiple IR, cm, and mm sources are found at the flow interface as well (Figs. 1 and 4 of Juárez et al. 2017). Overall, the similarities of NGC 6334 V and ISOSS J23053+5953 are remarkable, despite clearly being located in different parts of the Galaxy. Converging flows are also reported in DR21(OH) (Csengeri et al. 2011), Orion (Lee et al. 2013), G035.39-00.33 (Jiménez-Serra et al. 2014), G0.253+0.016 (Johnston et al. 2014), the Galactic Center (Schwörer et al. 2019), the Musca filament (Bonne et al. 2020), and G+0.693-0.03 (Zeng et al. 2020).

Toward core 3 in ISOSS J23053+5953, a chemical differentiation between DCO⁺ peaking toward the north (Fig. A.2) and N₂ D⁺ peaking toward the west (Fig. A.21) of the mm core is observed. In general, N₂ D⁺ is less depleted than DCO⁺ due to N₂ residing longer in the gas phase compared to CO and CO leading to an enhancement of deuterated species (Emprechtinger et al. 2009). Due to the fact that all CO isotopologues show wide spread emission, the difference between the two species might not be due to CO freeze-out. Murillo et al. (2018) also find spatially offset emission of these two species in the cold envelope of the low-mass protostar IRAS 16293-2422. In their models, the different locations of the molecular emission are due to a temperature effect, where DCO⁺ is located in a colder region at the disk-envelope interface. The colliding flow might have caused a temperature change toward core 3 or the spatial difference might be caused by potential molecular outflows. Two jet-like features are observed in SiO (Jet1 and Jet2), both pointing to core 3 (Fig. 7).

7 Conclusions

We used high angular resolution observations at 1.3 mm of the young star-forming regions ISOSS J22478+6357 and ISOSS J23053+5953 and MIR and FIR archival data from the Spitzer and Herschel telescopes to study the physical and chemical properties of fragmented mm cores. This study (CORE-extension) is a follow-up of the CORE large program, for which the fragmentation and physical-chemical properties of fragmented cores toward HMSFRs classified as HMPOs, HMCs, and UCHII regions are presented in Beuther et al. (2018) and Gieser et al. (2021), respectively. Here, we focused on the physical (temperature and density profiles, outflows) and chemical properties of younger mm cores.

Our main results are summarized as follows:

1. Comparing the 1.3 mm continuum with Spitzer IRAC and MIPS observations, we find that some of the fainter mm cores already have strong MIR counterparts suggesting clustered subsequent star formation occurring in both regions. A generation of more evolved low-mass protostars surrounds the younger intermediate- and high-mass protostars. The bolometric luminosity of the cores, which can be resolved in the MIR and FIR, was estimated from fitting their SEDs from MIR to mm wavelengths. Two or three temperature components were necessary to properly fit the full SEDs.

2. The radial density profile was inferred from the 1.3 mm continuum visibilities and the radial temperature profile was inferred from H₂CO rotation temperature maps. We assumed that both profiles follow a power-law profile, and we derive a mean density and temperature power-law index of p = 1.6 ± 0.4 and q = 0.4 ± 0.3, respectively. We find similar values toward the cores of the original more evolved CORE sample (Gieser et al. 2021), suggesting that both the radial density and temperature power-law index do not change significantly from the IRDC to HMC stage.

3. The 1.3 mm spectra are generally line-poor, compared to spectra observed in the original CORE sample, and mostly simple species with low upper energy levels are detected. Molecular column densities were inferred with XCLASS. The abundances (relative to C¹⁸O) of the cores compared to the cores of the original CORE sample (Gieser et al. 2021) are about one order of magnitude lower suggesting that most species are still depleted in the younger ISOSS sources studied here.

4. The CO (2−1) line wing emission reveals multiple bipolar molecular outflows in both regions. In ISOSS J22478+6357, we find five mm cores with outflows that are aligned almost parallel in the plane of the sky. These mm cores are connected by a filamentary structure seen in DCO⁺ and the outflows are perpendicular to this filament. A large-scale bipolar outflow is seen toward core 1 in ISOSS J23053+5953, which is a high-mass protostar.

5. Most of the mm cores do not show distinct molecular emission peaks (except for cores 1 and 2 in ISOSS J22478+6357, and for cores 1, 2, and 6 in ISOSS J23053+5953). Multiple molecular emission peaks are detected within the FOV in both regions with no associated mm continuum cores for transitions of SO, SiO, H₂CO, and CH₃OH being well-known shock tracers. These regions can be linked to bipolar outflows ISOSS J22478+6357 and a combination of a colliding flow and outflows in ISOSS J23053+5953 creating shocks in the gas surrounding the mm cores.

6. We applied the statistical tool HOG to the line data in both regions to find velocity channels with a high spatial correlation between two molecular transitions. As expected, transitions of the same species have a high correlation, but also the shock-tracing species and the molecules tracing the larger-scale emission are spatially correlated with no significant velocity offsets. With HOG it will be possible for future programs to easily identify spatial correlations in a large sample of regions and many detected transitions on a channel-by-channel comparison.

7. A large-scale velocity gradient is seen in ISOSS J23053+5953 toward the cores 2, 3 and 6. Between these cores, a collidingflow is creating a shock enhancing the abundances of many detected species. This colliding flow might have triggered further star formation.

8. With the physical-chemical model MUSCLE, we estimate chemical timescales of four cores with the highest number of detectedspecies and for which we were able to derive density and temperature profiles. We find a mean chemical timescale of 20 000 yr, while for the cores of the original more evolved CORE sample a mean of 60 000 yr is derived (Gieser et al. 2021). Most cores in the ISOSS regions are in an early evolutionary phase, which is supported by the shorter chemical timescales, smaller molecular abundances, and lower temperatures. However, we also find some more evolved cores in the region with higher temperatures, signposts of complex chemistry and MIR signatures that indicates subsequent on-going star formation (cores 1 and 6 in ISOSS J23053+5953).

9. Comparing the mm core properties studied in the ISOSS and CORE regions, we find that some properties clearly change with time. The 1.3 mm peak intensity and integrated flux, core mass, H₂ column density, and kinetic temperature of a mm core increase with time, since these properties are anticorrelated with the M/L ratio and correlated with the chemical timescale τ_chem.

The CORE and CORE-extension studies reveal that subarcsecond angular resolution observations are required in order to investigate the properties of fragmented cores that show complex dynamical and chemical properties on these spatial scales (<10 000 au). The number of observed and analyzed mm cores is still not high; however, this analysis demonstrates a pathway for ongoing and future large programs targeting hundreds or thousands of embedded young (high-mass) protostars. The 1.3 mm setup of the CORE project is not ideal for studying the molecular composition in young and cold regions in the IRDC phase, as well as in evolved UCHII regions. High angular resolution observations at 3 mm covering the ground-state transitions of many simple species, including many deuterated molecules, would provide a better insight on the earliest evolutionary stages. An intriguing result of this study is that both regions harboring young intermediate- and high-mass protostars are already embedded in a cluster of more evolved low-mass MIR-bright YSOs.

Acknowledgements

The authors would like to thank the anonymous referee whose comments helped to improve the clarity of this paper. 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). C.G. and H.B. acknowledge support from the European Research Council under the Horizon 2020 Framework Programme via the ERC Consolidator Grant CSF-648505. D.S. acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: “Building a Habitable Earth (SE 1962/6-1)”. R.K. acknowledges financial support via the Emmy Noether and Heisenberg Research Grants funded by the German Research Foundation (DFG) under grant no. KU 2849/3 and 2849/9. A.P. acknowledges financial support from the UNAM-PAPIIT IN111421 grant and the Sistema Nacional de Investigadores of CONACyT. R.G.-M. acknowledges support from UNAM-PAPIIT project IN104319. A. P. and R.G.-M. also acknowledge support from CONACyT Ciencia de Frontera project number 86372. R.E.P. is supported by an NSERC Discovery grant. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). This research made use of photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2021).

Appendix A Moment maps

In Sect. 4.2.1 we used the line integrated intensity in order to study and compare the spatial morphology of the molecular emission. The moment 0 map of H₂CO 3_0,3 − 2_0,2, SiO 5 − 4, and CH₃OH 4_2,3 − 3_1,2E are shown in Sect. 4.2.1 in Figs. 6, 7, and 8, respectively. Here, the line integrated intensity maps of the remaining transitions (Table 2) are presented in Figs. A.1 −A.23. The line emission was integrated from v_LSR − 3 km s⁻¹ to v_LSR + 3 km s⁻¹. The intensity-weighted peak velocity (moment 1) map of molecular lines with spatially extended emission in ISOSS J23053+5953 is shown in Fig. A.24.

Fig. A.1 Thesame as Fig. 6, but for SO 55 − 44.

Fig. A.2 Thesame as Fig. 6, but for DCO+ 3 − 2.

Fig. A.3 Thesame as Fig. 6, but for H2S 22,0 − 21,1.

Fig. A.4 Thesame as Fig. 6, but for CH3OH 51,4 − 42,3E.

Fig. A.5 Thesame as Fig. 6, but for c-C3H2 61,6 − 50,5.

Fig. A.6 Thesame as Fig. 6, but for HC3N 24 − 23.

Fig. A.7 Thesame as Fig. 6, but for H2CO 32,2 − 22,1.

Fig. A.8 Thesame as Fig. 6, but for H2CO 32,1 − 22,0.

Fig. A.9 Thesame as Fig. 6, but for OCS 18 − 17.

Fig. A.10 The same as Fig. 6, but for C18O 2 − 1.

Fig. A.11 The same as Fig. 6, but for HNCO 100,10 − 90,9.

Fig. A.12 The same as Fig. 6, but for H${}_2^{13}$CO 31,2 − 21,1.

Fig. A.13 The same as Fig. 6, but for SO 65 − 54.

Fig. A.14 The same as Fig. 6, but for CH3OH 80,8 − 71,6E.

Fig. A.15 The same as Fig. 6, but for 13CO 2 − 1.

Fig. A.16 The same as Fig. 6, but for CH3CN 120 − 110.

Fig. A.17 The same as Fig. 6, but for CH3OH 81,8 − 70,7E.

Fig. A.18 The same as Fig. 6, but for CO 2 − 1.

Fig. A.19 The same as Fig. 6, but for OCS 19 − 18.

Fig. A.20 The same as Fig. 6, but for 13CS 5 − 4.

Fig. A.21 The same as Fig. 6, but for N2D+ 3 − 2.

Fig. A.22 The same as Fig. 6, but for SO2 42,2 − 31,3.

Fig. A.23 The same as Fig. 6, but for HC3N 26 − 25.

Fig. A.24 Intensity-weighted peak velocity (moment 1) maps of molecular lines with extended emission in ISOSS J23053+5953. The moment 1 map is presented in color scale (using a threshold of S∕N ≥ 3 of the line integrated intensity). The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar.

Appendix B Molecular column densities derived with XCLASS

In Sect. 4.4 we present H₂CO parameter maps derived with XCLASS. The calculation of the modeled spectral lines in XCLASS follows the assumption of LTE and a 1D isothermal source. In the dense parts of star-forming regions the LTE conditions are valid. Using the myXCLASSMapFit function, we fitted in each pixel the observed H₂CO spectral line data with one emission component within the full FOV if the peak intensity is >1 K. This threshold corresponds to a minimum S/N of ~5 in the high-resolution data (Table 3). In both regions, toward the edges of the mosaic the noise increases and therefore the fit results toward the edges of the mosaic are not reliable.

In XCLASS, the best-fit parameters can be obtained using several algorithms with the implemented model optimizer package MAGIX (Möller et al. 2013). We used an algorithm chain utilizing the Genetic algorithm, optimizing global minima, and the Levenberg-Marquart algorithm, optimizing local minima, with 50 iterations each to derive the best-fit parameter set for each species. Thetransition properties of the fitted lines are taken from the Virtual Atomic and Molecular Data Centre (VAMDC, Endres et al. 2016) with entries from the Cologne Database for Molecular Spectroscopy (Müller et al. 2005, CDMS,) and the Jet Propulsion Laboratory (JPL, Pickett et al. 1998) catalog. The line properties of the three fitted H₂CO transitions are summarized in Table 2. We varied each fit parameter within the following ranges: θ_source = 0.1 − 2″, T_rot = 1 − 200 K, N = 10¹² − 10¹⁷ cm⁻², Δv = 0.5 − 10 km s⁻¹, and v_off = −10 − +10 km s⁻¹. A further input is the beam size of the interferometric observations θ_beam and as in XCLASS only a single value can be given as an input, we computed the mean of the major and minor axis of the synthesized beam: ${\theta }_{\text{beam}}=\frac{{\theta }_{\text{maj}}+{\theta }_{\text{min}}}{2}$ (the synthesized beam of each spectral line data product is listed in Table 3).

In order to quantify the molecular content of all mm cores (listed in Table 4), we fitted the spectral line data with XCLASS in order to derive molecular column densities (Sect. 4.6). As we analyze the spectrum at the position of the core and not the full FOV of both regions, minor changes were employed in the XCLASS fitting routine in order to make the fit results as reliable as possible: In order to determine a robust velocity v_LSR for each core, we first fitted the C¹⁸O 2 − 1 transition with one emission component. We then corrected all spectra with the obtained velocity offset. The results for the systemic velocity of each core are summarized in Table 4. The fit parameter ranges for the source size and rotation temperature were θ_source = 0.1 − 10″, T_rot = 1 − 200 K, respectively. Setting 1 K as the lower rotation temperature limit is in no opposition to the fact that the gas temperatures should be higher than the cosmic microwave background (CMB), T_CMB = 2.73 K, as toward all species and cores the derived rotation temperature is T_rot > 3 K. The fit parameter ranges for the column density N and line width Δv were specified individually for each molecule. The emission of each line was fitted around ν₀ ± 6 MHz (the rest frequency ν₀ of all transitions is listed in Table 2). With both regions being line-poor in general, most of the detected lines do not suffer from line blending, except for CH₃CN. The transition of CH₃CN 12₀ − 11₀ is blended with 12₁ − 11₁. Both transitions were considered when fitting the molecular emission in XCLASS.

The same algorithm chain as used for the H₂CO parameter maps for the was used. In order to estimate the uncertainties of the fit parameters, the Markov Chain Monte Carlo (MCMC) error estimation algorithm was used afterward with 50 iterations.

The column densities of all fitted species and 29 cores are listed in Tables B.1 and B.2. The uncertainties were estimated using the MCMC error estimation algorithm in XCLASS. Upper limits are listed for fits that were discarded. The constraints are explained in Sect. 4.6. An example of the observed and XCLASS modeled spectrum for all fitted lines of core 1 in ISOSS J22478+6357 and ISOSS J23053+5953 is shown in Fig. B.1 and Fig. B.2, respectively.

Fig. B.1 Example spectrum of all molecular transitions fitted with XCLASS toward core 1 in ISOSS J22478+6357. In each panel, the observed spectrum of the emission line and the corresponding XCLASS fit are presented by the black and red lines, respectively. The green dotted line shows the 3σline level of the transition (Table 3). The gray dashed line indicates the rest frequency of the transition (Table 2). If XCLASS fits were discarded no fit is shown (the constraints are explained in Sect. 4.6). An example spectrum of core 1 in ISOSS J23053+5953 is shown in Fig. B.2.

Fig. B.2 The same as Fig. B.1, but for core 1 in ISOSS J23053+5953.

Appendix C MUSCLE results

For four cores in our sample, which have enough species detected in the 1.3 mm setup and for which the density and temperature profiles could be derived, we applied the physical-chemical model MUSCLE in Sect. 5 to estimate the chemical timescale τ_chem. As the initial chemical composition of the gas is unknown and the physical structure is static in MUSCLE, we modeled each core with three initial chemical conditions, referred to as the HMPO, HMC, and UCHII model. The initial chemical conditions are taken from a study of 59 HMSFRs by Gerner et al. (2014, 2015).

In each time step, the modeled column densities, computed with ALCHEMIC, are compared with the observed column densities. The time step with the lowest χ² is taken as the best fit. Table C.1 shows for each of the four modeled cores and each initial condition model, the best-fit chemical timescale τ_chem, the corresponding χ² value, and the percentage of well modeled species ϒ. A molecule is considered as “well modeled” when the modeled and observed column density agree within a factor of ten. A weighted ${\overline{\chi }}^2$ was computed taking into account the percentage of well modeled species ϒ: ${\overline{\chi }}^2=\frac{{\chi }^2}{{{\rm Y}}^2}$. For each core, the initial condition model with the lowest ${\overline{\chi }}^2$ was taken as the best-fit model.

References

All Tables

All Figures

Fig. 1 NOEMA 1.3 mm continuum emission of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The NOEMA continuum data are presented in color scale and pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. The core outer radius rout is indicated by a dashed blue circle. The synthesized beam is shown in the bottom left corner. A linear spatial scale of 0.1 pc is indicated by a black scale bar.

In the text

Fig. 2 Multi-wavelength overview of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The top panels show a Spitzer IRAC composite RGB image at 4.5 (blue), 5.8 (green), and 8.0 μm (red). Emission with S∕N < 5 above the sky background is clipped and shown in black. The middle panels show Spitzer MIPS data at 24 μm. The bottompanels show a composite RGB image of Herschel PACS observations at 70 (blue), 100 (green), and 160 μm (red). Emission with S∕N < 5 above the sky background is clipped and shown in black. The Herschel clumps studied in Ragan et al. (2012) are indicated by green circles and labeled in green. In all panels, the NOEMA continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The area outside of the NOEMA field-of-view is gray-shaded. All mm cores identified in Beuther et al. (2021) are labeled in pink. The synthesized beam of the NOEMA continuum data is shown in the bottom left corner. A linear spatial scaleof 0.1 pc is indicated by a white scale bar.

In the text

Fig. 3 Radial visibility profiles of cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953. The azimuthally averaged visibility amplitudes $\overline{V}$, considering the real and imaginary components, are shown in black. The bottom x-axis shows the uv distance and the top x-axis is converted to the corresponding linear spatial scale. A fit to the data and its uncertainties (±1σ) are shown by the red solid and dashed lines, respectively.

In the text

	Fig. 4 SED of mm cores with associated Spitzer IRAC emission in at least three bands. The fluxes (black points) are obtained from Spitzer IRAC (3.6, 4.5, 5.8, and 8.0 μm), Spitzer MIPS (24 μm), Herschel PACS (70, 100, and 160 μm), and NOEMA 1.3 mm observations (Tables 4 and 6). The total SED fit is shown by a red line and individual components are indicated by blue lines.
In the text

Fig. 5 CO 2−1 outflows in ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The blue and red contours show the integrated intensity of the blue- and red-shifted line wings, respectively. The intensity is integrated over vLSR − 21 to vLSR − 6 and over vLSR + 6 to vLSR + 21 km s−1 for the blue- and red-shifted integrated intensity, respectively. The corresponding velocity ranges are shown in each panel. Contour levels range from 10 to 90% of the peak integrated intensity with steps of 10%. In ISOSS J22478+6357 the red- and blue-shifted peak integrated intensity is 45 and 47 K km s−1 and in ISOSS J23053+5953 it is 112 and 140 K km s−1, respectively. Red- and blue-shifted outflow directions are indicated by red and blue arrows, respectively. The NOEMA 1.3 mm continuum data are shown in black contours. The dotted black contour marks the − 5σcont level. Solid black contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in black. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). Jet-like structures, seen in SiO 5− 4 emission, are indicated by orange ellipses (Fig. 7). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated bya black scale bar.

In the text

Fig. 6 Line integrated intensity map of H2CO 30,3 − 20,2 of ISOSS J22478+6357 (left panel) and ISOSS J23053+5953 (right panel). The integrated intensity of the transition with a threshold of S∕N ≥ 3 is presented in color scale. The line intensity is integrated from a velocity of vLSR − 3 km s−1 to vLSR + 3 km s−1. It is marked if the transition is detected (√) or not (×), further explained in Sect. 4.2.1. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). Jet-like structures, seen in SiO 5− 4 emission, are indicated by orange ellipses (Fig. 7). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated bya black scale bar.

In the text

	Fig. 7 Same as Fig. 6, but for SiO 5−4.
In the text

	Fig. 8 Same as Fig. 6, but for CH3OH 42,3 − 31,2E.
In the text

	Fig. 9 HOG results for ISOSS J22478+6357 (top panel) and ISOSS J23053+5953 (bottom panel). The peak projected Rayleigh statistic V is shown for each transition pair (Sect. 4.2.2). Examples of high spatial correlations in velocity space are shown in Fig. 10.
In the text

Fig. 10 HOG correlation results for DCO+ (3− 2) and H2CO (30,3 − 20,2) toward ISOSS J22478+6357 (top panel) and for SiO (5−4) and CH3OH (42,3 − 31,2E) toward ISOSS J23053+5953 (bottom panel). Emission maps of the channels with the highest correlation are shown in the left and middle panels. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated by a black scale bar. The projected Rayleigh statistic V, computed for all velocity pairs, is shown in the right panel. The red dashed line indicates equal velocity in both transitions. For ISOSS J23053+5953, the velocity regime around the two velocity components seen in line emission at v ≈ −52 km s−1 and − 50 km s−1 (Sect. 4.3) are highlighted by black dash-dotted circles.

In the text

Fig. 11 H2CO parameter maps of ISOSS J22478+6357 derived with XCLASS. The fit parameter maps (Trot, N, Δv, voff) for H2CO are presented in color scale. The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. All mm cores identified in Beuther et al. (2021) are labeled in pink. Positions with a peak in molecular emission, but no associated continuum (“molecular peaks”, MP) are indicated by orange circles (Sect. 4.2.1). The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar. A linear spatial scale of 0.1 pc is indicated by a black scale bar.

In the text

	Fig. 12 Same as Fig. 11, but for ISOSS J23053+5953. Jet-like structures, seen in SiO 5−4 emission, are indicated by orange ellipses (Fig. 7).
In the text

Fig. 13 Radial temperature profile for cores 1−4 in ISOSS J22478+6357 and cores 1−6 in ISOSS J23053+5953. The radial H2CO temperature profile is shown by the black data points. The outer radius rout, estimatedwith the clumpfind algorithm in the continuum data (Table 4), is indicated by the dash-dotted purple vertical line. A fit to the data and its uncertainties (±1σ) are shown by the red solid and dashed lines, respectively.

In the text

	Fig. 14 Abundance histograms (relative to C18O) of all species fitted with XCLASS. The gray histogram shows the combined results of all 29 cores. The mean and standard deviation are shown in the top left corner. The green dashed and blue dash-dotted histograms show the results for the cores in ISOSS J22478+6357 and ISOSS J23053+5953, respectively.
In the text

	Fig. 15 Comparison of the observed and modeled column densities of all cores modeled with MUSCLE. The observed column density is shown in red and upper limits are indicated by arrows (Tables B.1 and B.2). The column densities of the MUSCLE best-fit model are shown by the blue data points assuming a factor of ten uncertainty.
In the text

	Fig. 16 Mean abundances (relative to C18O) and standard deviation of the cores in both ISOSS regions (green circles) in comparison with the cores in the CORE sample (blue triangles) taken from Gieser et al. (2021). Only species that were studied both in this work and in Gieser et al. (2021) are shown.
In the text

Fig. 17 Core parameters with a Spearman correlation coefficient of |rS|≥ 0.4. The results for the cores in ISOSS J22478+6357 and ISOSS J23053+5953 are shown in green (this work, Tables 4–6) and the results for the cores of the original CORE sample are presented in blue (Gieser et al. 2021). In each panel, the correlation coefficient rS and p-value are shown in the top left corner. For completeness, the clump-average M/L ratio taken from Ragan et al. (2012) of the corresponding ISOSS cores (Table 4) is shown in red.

In the text

	Fig. A.1 Thesame as Fig. 6, but for SO 55 − 44.
In the text

	Fig. A.2 Thesame as Fig. 6, but for DCO+ 3 − 2.
In the text

	Fig. A.3 Thesame as Fig. 6, but for H2S 22,0 − 21,1.
In the text

	Fig. A.4 Thesame as Fig. 6, but for CH3OH 51,4 − 42,3E.
In the text

	Fig. A.5 Thesame as Fig. 6, but for c-C3H2 61,6 − 50,5.
In the text

	Fig. A.6 Thesame as Fig. 6, but for HC3N 24 − 23.
In the text

	Fig. A.7 Thesame as Fig. 6, but for H2CO 32,2 − 22,1.
In the text

	Fig. A.8 Thesame as Fig. 6, but for H2CO 32,1 − 22,0.
In the text

	Fig. A.9 Thesame as Fig. 6, but for OCS 18 − 17.
In the text

	Fig. A.10 The same as Fig. 6, but for C18O 2 − 1.
In the text

	Fig. A.11 The same as Fig. 6, but for HNCO 100,10 − 90,9.
In the text

	Fig. A.12 The same as Fig. 6, but for H${}_2^{13}$CO 31,2 − 21,1.
In the text

	Fig. A.13 The same as Fig. 6, but for SO 65 − 54.
In the text

	Fig. A.14 The same as Fig. 6, but for CH3OH 80,8 − 71,6E.
In the text

	Fig. A.15 The same as Fig. 6, but for 13CO 2 − 1.
In the text

	Fig. A.16 The same as Fig. 6, but for CH3CN 120 − 110.
In the text

	Fig. A.17 The same as Fig. 6, but for CH3OH 81,8 − 70,7E.
In the text

	Fig. A.18 The same as Fig. 6, but for CO 2 − 1.
In the text

	Fig. A.19 The same as Fig. 6, but for OCS 19 − 18.
In the text

	Fig. A.20 The same as Fig. 6, but for 13CS 5 − 4.
In the text

	Fig. A.21 The same as Fig. 6, but for N2D+ 3 − 2.
In the text

	Fig. A.22 The same as Fig. 6, but for SO2 42,2 − 31,3.
In the text

	Fig. A.23 The same as Fig. 6, but for HC3N 26 − 25.
In the text

Fig. A.24 Intensity-weighted peak velocity (moment 1) maps of molecular lines with extended emission in ISOSS J23053+5953. The moment 1 map is presented in color scale (using a threshold of S∕N ≥ 3 of the line integrated intensity). The NOEMA 1.3 mm continuum data are shown in pink contours. The dotted pink contour marks the − 5σcont level. Solid pink contour levels are 5, 10, 20, 40, and 80σcont. The synthesized beam of the continuum data is shown in the bottom left corner. The synthesized beam of the spectral line data is similar.

In the text

Fig. B.1 Example spectrum of all molecular transitions fitted with XCLASS toward core 1 in ISOSS J22478+6357. In each panel, the observed spectrum of the emission line and the corresponding XCLASS fit are presented by the black and red lines, respectively. The green dotted line shows the 3σline level of the transition (Table 3). The gray dashed line indicates the rest frequency of the transition (Table 2). If XCLASS fits were discarded no fit is shown (the constraints are explained in Sect. 4.6). An example spectrum of core 1 in ISOSS J23053+5953 is shown in Fig. B.2.

In the text

	Fig. B.2 The same as Fig. B.1, but for core 1 in ISOSS J23053+5953.
In the text