Physical and chemical structure of high-mass star-forming regions

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Physical and chemical structure of high-mass star-forming regions C. Gieser, H. Beuther, D. Semenov, A. Ahmadi, S. Suri, T. Möller, M. T Beltrán, P. Klaassen, Q. Zhang, J. S Urquhart, et al.


Introduction
The development of large telescopes and highly sensitive instruments has provided the opportunity to investigate star-forming regions within and even outside of the Milky Way in great detail (e.g., Larson 1981;Shu et al. 1987;Kennicutt 1998;McKee & Ostriker 2007). The study of high-mass star formation (HMSF) is challenging from an observational point of view (e.g., Beuther et al. 2007a;Bonnell 2007;Zinnecker & Yorke 2007;Smith et al. 2009;Tan et al. 2014;Krumholz 2015;Schilke 2015;Motte et al. 2018;Rosen et al. 2020). High-mass stars are less commonly observed as they are typically located at large distances of several kiloparsecs and therefore difficult to spot at a high spatial  are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5)  resolution. The evolution of HMSF is fast, that is, on the order of 10 5 yr, as revealed by observations (see, e.g., Table 2 in Motte et al. 2018;Mottram et al. 2011a) and theoretical models (McKee & Tan 2002Kuiper & Hosokawa 2018). Therefore, highmass protostars are still deeply embedded within their parental molecular cloud when they reach the main sequence.
In contrast to low-mass star formation, no clear consensus has yet been reached on the evolutionary stages of high-mass stars (see, e.g., Fig. 14 in Motte et al. 2018). High-mass starforming regions (HMSFRs) can be roughly categorized into several evolutionary stages based on their observed properties at infrared and radio wavelengths. Various classifications exist in the literature, based on, for instance, infrared properties. Cooper (2013) classifies massive young stellar objects (MYSOs) into three types: type I showing strong H 2 , but no ionized lines; type II being weaker in H 2 , but showing HI emission in the Brackett series; type III having strong HI and weak H 2 line emission. These near-infrared observations show that type III MYSOs are bluer than type I MYSOs, indicating that they are more evolved. A&A 648, A66 (2021) In this study, we follow the observationally driven nomenclature of four different evolutionary stages during HMSF by Gerner et al. (2014Gerner et al. ( , 2015. When regions of a molecular cloud undergo gravitational collapse and fragmentation, protostars form in cold and dense clumps, often detected as infrared dark clouds (IRDCs, e.g., Rathborne et al. 2006). These clouds have typical temperatures of ∼10−20 K and are visible at (sub)mm wavelengths (e.g., Pillai et al. 2006), but show none or only weak emission in the infrared. As the collapse continues, due to the conservation of angular momentum, a disk-outflow system forms around the protostar and the temperature of the surrounding envelope increases due to central heating by the protostar (e.g., Sánchez-Monge et al. 2013a;Beltrán & de Wit 2016). Through the accretion of the infalling gas, protostars increase their masses and luminosities, becoming visible at infrared wavelengths. Protostars during this stage are referred to as high-mass protostellar objects (HMPOs, e.g., Sridharan et al. 2002;Williams et al. 2004;Beuther et al. 2010;Duarte-Cabral et al. 2013). Such HMPOs are characterized by their high bolometric luminosities, L > 10 4 L , a strong thermal dust continuum, but only weak cm emission (Sridharan et al. 2002;Beuther et al. 2002). As the temperature of the surrounding envelope reaches 100 K, molecules that had formed on the surfaces of dust grains and resided in ice layers evaporate into the gas-phase revealing a rich molecular reservoir (e.g., Osorio et al. 1999). In this stage, the objects are referred to as hot molecular cores (HMCs) or "hot cores" (e.g., Cesaroni et al. 1997). The HMPOs and HMCs have similar physical properties, however, the chemical composition of the gas around HMCs is richer due to higher temperatures in the envelope. Low-mass analogues are referred to as "hot corinos" (e.g., IRAS 16293-2422, Bottinelli et al. 2004). As the mass growth continues, the strong protostellar radiation dissociates and ionizes the surrounding envelope gas and a hyper or ultra-compact HII (HC/UCHII) region forms (e.g., Wood & Churchwell 1989;Hatchell et al. 1998;Churchwell 2002;Sánchez-Monge et al. 2013b). HII regions can be studied and classified through their strong free-free emission at cm wavelengths (e.g., Peters et al. 2010). These four evolutionary stages are not sharp transitions and overlap, especially between the intermediate HMPO/HMC stages, and there are HMCs that have already formed a HC/UCHII region.
Aside from the physical complexity of HMSFRs, observations reveal one of the most complex chemical compositions in the interstellar medium (ISM, for a review, see, e.g., Herbst & van Dishoeck 2009;Jørgensen et al. 2020). To date, about 200 molecules have been detected in the ISM (an overview is found in McGuire 2018) and most of these can be observed toward the high-mass star-forming molecular cloud Sagittarius B2 in the Galactic center (e.g., Belloche et al. 2013). Spectroscopic studies at mm and sub-mm wavelengths have revealed the composition of the molecular gas in HMSFRs to be extremely diverse: nitrogen (N)-bearing species (such as HCN, CH 3 CN); oxygen (O)-bearing species (such as H 2 CO); sulfur (S)-bearing species (such as SO, SO 2 ). In shocked regions, enhanced abundances of S-bearing molecules such as SO 2 and silicon (Si)-bearing molecules such as SiO have been observed (Schilke et al. 1997). During the HMC stage, a large variety of so-called complex organic molecules (COMs) evaporate into the gas-phase and produce line-rich spectra at mm wavelengths. We use the definition by Herbst & van Dishoeck (2009) stating that a COM contains six or more atoms. In the densest regions where protostars form, the chemical composition of the molecular gas has been studied thoroughly. An important finding is the chemical segregation of O-and N-bearing species (e.g., Rodgers & Charnley 2001). This has been observed and studied toward many bright HMCs, such as Orion-KL (Caselli et al. 1993;Feng et al. 2015); W3 H2O and W3 OH (Wyrowski et al. 1999;Qin et al. 2015); AFGL 2591 (Jiménez-Serra et al. 2012;Gieser et al. 2019); NGC 7538 IRS9, W3 IRS5 and AFGL 490 (Fayolle et al. 2015); G35.20-0.74N (Allen et al. 2017); SgrB2(N) (Bonfand et al. 2017;Mills et al. 2018); or AFGL 4176 .
In this paper, we want to pin down the physical and chemical properties of HMSFRs from the dust and molecular line emission and use this information to investigate the chemical timescales. We analyze the physical structure and molecular content of 18 HMSFRs in combination with a physical-chemical model on spatial scales ranging from 10 000 au down to our resolution limit of 300 au. We successfully tested the method presented in this paper on the well-studied hot core AFGL 2591 VLA3 (Gieser et al. 2019). In this case study, we used the molecular line and dust continuum data from the observations of the CORE project to derive the density and temperature structure and molecular column densities of this hot core and we applied a physical-chemical model to estimate the chemical age. Here, we apply this method to all 18 CORE regions with the goal of comparing and understanding the chemical diversity we observe during HMSF.
The paper is organized as followed: in Sect. 2, we summarize the observations, data calibration, and imaging techniques. In Sect. 3, we derive the temperature and density structure of the regions. The molecular content is analyzed in Sect. 4. We apply a physical-chemical model to the observed molecular column densities in Sect. 5. The results of the physical and chemical properties of the regions are discussed in Sect. 6. Our summary and conclusions are given in Sect. 7.

Observations
The NOrthern Extended Millimeter Array (NOEMA) large program known as "Fragmentation and disk formation during highmass star formation -CORE" is a high angular-resolution survey (∼0. 4 at 1.37 mm) designed to study the fragmentation, kinematic, and chemical properties of a homogeneous sample of 18 HMSFRs. The observations consist of spectral line and continuum data in Band 3 (1 mm). An overview of the project and results of the analysis of the dust continuum and fragmentation properties are presented in Beuther et al. (2018).
In combination with this paper, we release all data products of the CORE project with and without self-calibration conducted with GILDAS. The data can be found on the CORE collaboration website 1 . The CLEANed data products as well as the calibrated uv-tables are available at this link.

Sample
The sample of the CORE project consists of well-studied HMS-FRs on the northern hemisphere. They were selected based on their bolometric luminosity (L > 10 4 L ) and distance (d < 6 kpc), allowing us to study the evolution at early stages of HMSF in the HMPO-HMC stage. An overview of the properties of the regions is shown in Table 1 in this paper and additional properties, such as the luminosity L and mass M, are listed in Table 1 in Beuther et al. (2018). With declinations higher than +20 • , most of the regions are difficult or impossible to observe with the Atacama Large Millimeter/submillimeter Array (ALMA) and, therefore, can only be studied with Notes. The isotopic ratios are taken from Wilson & Rood (1994), further explained in Sects. 3.2, 4.2, and 5.1.
NOEMA at a high angular resolution (<1 ) at mm wavelengths. The angular resolution for all regions is homogeneous (∼0. 4), but the resolved linear spatial scales vary from 300 to 2300 au, as the regions are located at distances between 0.7 and 5.5 kpc. The regions show a large diversity of fragmentation properties : while some regions contain mainly a single isolated core (e.g., AFGL 2591), other regions fragment into up to 20 cores (e.g., IRAS 21078). Individual case-studies were carried out on some of the regions, for instance, the kinematic properties of W3 H2O , IRAS 23385 (Cesaroni et al. 2019), IRAS 23033 (Bosco et al. 2019), W3 IRS4 (Mottram et al. 2020), and IRAS 21078 (Moscadelli et al. 2021); the chemical composition of the pilot regions NGC 7538 S and NGC 7538 IRS1 (Feng et al. 2016), and AFGL 2591 (Gieser et al. 2019). A multi-wavelength modeling study of AFGL 2591 is presented in Olguin et al. (2020). In Ahmadi et al. (in prep.) the kinematic analysis of the complete CORE sample will be covered, while in this study, we focus on the physical structure (temperature and density) and chemical composition of the molecular gas. We do not include the pilot regions NGC 7538 S and NGC 7538 IRS1 presented in Beuther et al. (2018) in our analysis, as it is only for the remaining 18 regions, we have a homogeneous multi-configuration NOEMA data set. The pilot studies have no corresponding D array observations, and thus, less uv-coverage, so these could not be accurately compared. However, a detailed individual chemical analysis of the pilot regions is already presented in Feng et al. (2016).

NOEMA observations
The sample was observed from 2014 to 2017 in the A (extended), the old B or new C (intermediate), and D (compact) array configurations. Two regions were observed at a time in track-sharing pairs to reduce calibration time. Spectral line data were obtained with the broadband WideX correlator at a rest-frequency range of 217.2−220.8 GHz, with a spectral resolution of ∼2.7 km s −1 . In addition, eight narrow-band units were placed within this frequency range in order to obtain high spectral resolution data (∼0.4 km s −1 ) of kinematically interesting lines such as CH 3 CN (a summary of the high spectral resolution units is summarized in Tables 2 and 3 in Ahmadi et al. 2018). The interferometric data were calibrated with the CLIC package in GILDAS 2 . The 1.37 mm continuum data were extracted from the broadband WideX data by carefully selecting line-free channels in each region.

IRAM 30 m observations
Interferometers spatially filter the extended emission. The shortest baseline of the NOEMA array is ∼20 m, thus spatial scales larger than 16 are not recovered by the interferometric observations at 1.37 mm. All CORE regions were therefore observed with the IRAM 30 m telescope using the Eight MIxer Receiver (EMIR, Carter et al. 2012) in order to recover the missing flux in the spectral line data due to spatial filtering. EMIR has four basebands and each have a width of ∼4 GHz and a spectral resolution of 200 kHz corresponding to ∼0.3 km s −1 at 1.37 mm). The lower inner baseband (LI) of our EMIR spectral setup covers the same spectral range as the NOEMA observations with the broadband WideX correlator. The half power beam width (HPBW) of the IRAM 30 m telescope is 11. 8 at 1.37 mm. For a detailed description of the IRAM 30 m data calibration we refer to Appendix A in Mottram et al. (2020). The continuum data have no complementary single-dish observations, so here spatial filtering can remain an issue as discussed in Beuther et al. (2018).

Self-calibration
The majority of the continuum data have a high signal-to-noise ratio (S/N) in the GILDAS standard calibrated data (see Table 2).

Region
Synthesized beam GILDAS CASA GILDAS Continuum data Merged line data Standard calibration Self-calibration Self-calibration Notes. The noise of the continuum data σ cont is computed in a 3. 5 × 3. 5 rectangle in an emission-free region. I peak is the continuum peak intensity. The average map noise of the merged (combined NOEMA + IRAM 30 m) spectral line data σ line,map is computed in a line-free range from a rest-frequency of 219.00-219.13 GHz with a channel width of 3 km s −1 (2.2 MHz) within the full primary beam.
However, the phase noise can be high due to an unstable atmosphere during the observations that causes the flux to be scattered around the source. This issue can be improved by applying phase self-calibration to the interferometric data (e.g., Pearson & Readhead 1984;Radcliffe et al. 2016). The CORE continuum data presented in Beuther et al. (2018) were successfully phase self-calibrated using the Common Astronomy Software Applications package (CASA) through an iterative masking of the source. However, at that time, it was not possible to apply the self-calibration solution to the spectral line data. The CASA phase self-calibration of the continuum data was performed by applying an interactive mask in each self-calibration loop starting with the strongest structures first and proceeding with the weaker structures . Depending on the S/N of the region, we used solution intervals of 220, 100, or 45 s.
We used the self-calibration tool in GILDAS to phase selfcalibrate the continuum as well as the spectral line data of the CORE project to provide homogeneously calibrated data products of the full CORE observations. The crucial point of phase self-calibration is to start with a good enough spatial model of the source. The selfcal procedure of the GILDAS/MAPPING package uses as a source model the first CLEAN components n CLEAN found during a previous step of deconvolution. The basic idea is that the first CLEAN components deliver a model of the source with a high S/N whose spatial structure is not much affected by flux scattering. The number of CLEAN components n CLEAN must be large enough to get a fair representation of the source structure and small enough to avoid deconvolving scattered flux that would be confused with noise. The visibilities are usually averaged in time to increase their S/N during the first self-calibration iteration and this averaging time is progressively lowered to the minimum possible integration time in the following iterations. In practice, we started by deconvolving the continuum source with an absolute flux stopping criterion set to three times the continuum noise. This gives a number of CLEAN iterations, n iter , which we used to iterate the self-calibration three times, increasing n CLEAN from n iter /4, to n iter /2, and to n iter , and decreasing the averaging time from 200 s, to 100 s, and to 45 s. Only visibilities with a S/N > 6 were self-calibrated, but the remaining visibilities were kept in the proceeding CLEANing of the data in order to avoid losing visibilities of the longest baselines, which would decrease the angular resolution. With this method, simple structures, as well as regions with a complicated morphology, can all be successfully phase self-calibrated. As the continuum data has the highest sensitivity, we used the solution from the continuum self-calibration and applied the gain solution to the broadband and narrow-band spectral line data using the UV_CAL task.
An example of how self-calibration is increasing the quality of the continuum data quality is shown in Fig. 1 for W3 IRS4, which has a complex morphology in the continuum data. While the GILDAS standard calibrated data already reveals the complex structure of the region, the noise is high throughout the primary beam with many negative features. Applying the CASA and GILDAS self-calibration significantly lowers the noise and increases the peak intensity. The continuum noise and peak intensity are summarized for all regions in Table 2. The mean continuum noise, σ cont , is 0.74 mJy beam −1 in the standard calibrated data, while for the CASA and GILDAS self-calibrated data, the mean noise, σ cont , is lowered by ∼25% to 0.53 mJy beam −1 and 0.56 mJy beam −1 , respectively. In general, σ cont , is higher toward the strong continuum sources S106, CepA HW2, and W3 H2O. The mean S/N of the GILDAS standard calibrated continuum data is 74, while for the CASA and GILDAS selfcalibrated continuum data, the mean S/N is improved by a factor of two to 130 and 132, respectively. Even though both CASA and GILDAS self-calibration methods improve the data quality, due to the different techniques to define a source model (interactive masking and defining number of clean components, respectively), there are differences in the resulting images. In Fig. 1, a faint structure toward the NW of the continuum peak is seen in the GILDAS standard calibrated and self-calibrated image, but it is not recovered in the CASA self-calibrated image. Comparing the self-calibration results (Table 2), the peak intensities are higher in the GILDAS self-calibrated data, while a lower noise is achieved in the CASA self-calibrated data. For very faint structures and a complex source morphology, a careful interpretation of the self-calibrated data is recommended, but overall all the main features are recovered with both self-calibration methods.
A comparison between the standard and self-calibrated broadband spectral line data is shown in Fig. 2 for the CH 3 CN 12 3 −11 3 transition around the location of the 1.37 mm continuum peak in the W3 IRS4 region. The emission is less fuzzy and more compact in the self-calibrated data product. On individual spectra, self-calibration provides a significant increase of the line intensity which has a big impact, for example, when deriving column densities.

Data products and public release
The data merging of the interferometric and single-dish data and imaging is done in the GILDAS/MAPPING package. The WideX and the corresponding EMIR spectral-line data are smoothed to a common spectral resolution of 3.0 km s −1 . The narrow-band and the corresponding EMIR spectral line data are smoothed to a common spectral resolution of 0.5 km s −1 . In order to merge the NOEMA and IRAM 30 m data, the task UVSHORT converts the short spacings into a pseudo uv-table and combines them with the NOEMA data.
The deconvolution of the NOEMA continuum, NOEMA spectral line and merged (combined NOEMA + IRAM 30 m) spectral line data is performed in GILDAS/MAPPING using the Clark CLEAN algorithm (Clark 1980) and adopting three different weightings: robust weighting with a robust parameter of 0.1 (θ ∼ 0. 4); a robust parameter of 1.0 (θ ∼ 0. 6); and natural weighting (θ ∼ 1. 0). The stopping criterion for the continuum and broadband spectral line data is set to f res = 0.01, which Comparison of the GILDAS standard calibration (left) and GILDAS self-calibration (right) of the W3 IRS4 broadband spectral line data zoomed in toward the position of the continuum peak, shown in the top panels. The integrated intensity of the CH 3 CN 12 3 −11 3 transition is shown in color. In each panel, the beam size is shown in the bottom left corner and pink bar in the top left corner indicates a linear spatial scale of 5000 au. Bottom panel: CH 3 CN 12 3 −11 3 spectrum at the position of the W3 IRS4 continuum peak of the GILDAS standard calibrated data (black) and of the GILDAS self-calibrated data (green). The dashed orange line shows the systemic velocity of the region LSR and the dash-dotted blue lines show the lower and upper velocity limits ( LSR ± 10 km s −1 ) used for the integrated intensity map shown in the top panel.
corresponds to a minimum fraction of 1% of the peak intensity in the dirty image. The stopping criterion for the narrow-band spectral line data is either n iter = 5000 (maximum number of iterations) or a res = 0.01, which corresponds to a maximum intensity of 10 mJy beam −1 in the residual image.
A66, page 5 of 68 A&A 648, A66 (2021) Molecules such as CO isotopologues and SO may have extended emission in the merged data products (Mottram et al. 2020). As the Clark algorithm assumes emission from point sources, the CLEANed map may have point-like artifacts. The SDI algorithm (Steer et al. 1984) may improve the CLEANed image in such cases. A comparison between the Clark and SDI algorithms for imaging of the W3 IRS4 region is presented in Mottram et al. (2020). Final data products of the merged broadband spectral line data of molecular lines with potential large-scale emission ( 13 CO 2−1, SO 6 5 −5 4 , H 2 CO 3 0,3 −2 0,2 , H 2 CO 3 2,2 −2 2,1 , H 2 CO 3 2,1 −2 2,0 ) CLEANed using the SDI algorithm with robust weighting (robust parameter of 3, θ ∼ 0. 6), with a stopping criterion of n iter = 25 000 are provided as well. Primary beam correction is applied to the continuum and spectral line data products.
In this study, we use the primary beam corrected NOEMAonly continuum and merged (NOEMA + IRAM 30 m) broadband spectral line data. Both data products are imaged with the Clark algorithm and a robust parameter of 0.1, resulting in the highest angular resolution (θ ∼ 0. 4). Table 2 summarizes the properties of the standard and self-calibrated data products for all regions, including the synthesized beam (major axis θ maj , minor axis θ min , and position angle PA), noise of the continuum data σ cont , continuum peak intensity I peak , and noise in the merged (NOEMA + IRAM 30 m) spectral line data σ line,map . The mean map noise of the merged spectral line data is 0.44 K.

Physical structure
The linear spatial resolution of the CORE sample ranges from 300−2300 au. At this spatial resolution, it is not possible to resolve potential disks surrounding the protostars, but these can be studied and characterized through a kinematic analysis of the line profiles ). However, we do resolve the gas and dust envelopes around individual cores which can be approximated as spherically symmetric objects for which the radial temperature profile T (r) and radial density profile n(r) of the envelope gas can be described by power-laws (e.g., van der Tak et al. 2000;Beuther et al. 2002;Palau et al. 2014): with T in = T (r in ) and n in = n(r in ) at an arbitrary radius r in . The temperature power-law index q and density power-law index p of the core envelopes are important properties of HMS-FRs. By studying these quantities with observations, theoretical analytical models on how massive stars form can be constrained.

Continuum emission
A detailed analysis of the continuum data is presented in Beuther et al. (2018). The updated GILDAS self-calibrated continuum data are shown in contours in Fig. 3 In this study, we aim to analyze the physical properties and chemical variation across the regions, therefore, we selected a number of positions throughout the different regions. Deriving the column density of all detected species in every single spectrum in each region is computationally expensive and we restrict this method to the analysis of the H 2 CO and CH 3 CN temperature maps (Sect. 3.2). In order to study the cores and envelopes around forming protostars, we selected positions which have a clear spherically symmetric core-like morphology in the 1.37 mm continuum data which is the case for most regions toward the continuum peak position, but also multiple spherically symmetric objects within a region can be identified (e.g., toward IRAS 23033). In addition, we select positions in the extended envelopes to study potential differences in the chemical abundances. The in total, 120 selected positions are summarized in Table A.1 and marked in Fig. 3. The nomenclature of each position is denoted by the region name and increasing number with decreasing 1.37 mm continuum intensity. For each region, position 1 is toward the 1.37 mm continuum emission peak.
In order to estimate the H 2 column density (Sect. 4.1), we require for all positions to be detected in the 1.37 mm continuum (I 1.37 mm ≥ 5σ cont ). The number of selected positions is higher toward regions with extended envelopes and complex morphologies (e.g., IRAS 21078 and G138.2957) compared to compact regions (e.g., G139.9091 and S106). In this study, we define (in Sect. 3.2) a "core" as an object that shows a radially decreasing temperature profile over at least the width of two beams. In Sect. 5, we apply a 1D physical-chemical model to these defined cores using the observed temperature and density structure analyzed in Sects. 3.2 and 3.3, respectively, and observed molecular column densities (Sect. 4). These positions are labeled as "C" (core) in Table A.1. The remaining positions are locations in the dust envelope and environment around the cores. There are a few positions with a clear spherically symmetric core morphology in the dust emission, but no temperature profile can be derived (e.g., position 1 in S106). Cores with unresolved radial temperature profiles are also not included in this approach. In total we select 120 positions including 22 core positions. The broadband spectra of the 120 positions are shown in Fig. F.1. Table A.1 lists the noise of the broadband spectrum extracted from these positions and the systemic velocity LSR determined from the C 18 O 2−1 line, which may differ from the average region LSR (listed in Table 1) due to velocity gradients within the region (see Sect. 4.2). In contrast to this core definition, in the analysis by Beuther et al. (2018), "cores" are defined as fragmented objects with emission ≥10σ cont using the clumpfind algorithm (Williams et al. 1994) and, thus, more cores are found in their analysis.

Temperature structure
To reliably determine the rotation temperature, T rot , it is required to observe multiple transitions of a molecule. We use formaldehyde (H 2 CO) and methyl cyanide (CH 3 CN) as thermometers to determine the temperature structure of the regions, as both have multiple strong and optically thin lines in our spectral setup. The spectral line properties of these molecules are listed in Table B.1.
We model the spectral line emission of these molecules using the eXtended CASA Line Analysis Software Suite (XCLASS,  Table A.1) are marked and labeled in red (core positions) and blue (all remaining positions). The dashed black contours show the −5σ cont emission and the solid black contours start at 5σ cont with contour steps increasing by a factor of 2 (e.g., −5, 5, 10, 20, 40, 80,...σ cont ; see Table 2 for values of σ cont for each region). In each panel, the beam size is shown in the bottom left corner and the pink bar in the top left corner indicates a linear spatial scale of 5000 au. The field of view in each region is adjusted to only show the area with significant emission (≥5σ cont ), and for regions with wide-spread emission, the primary beam is indicated by a black circle. The broadband spectral setup includes three strong formaldehyde lines, two of which have the same upper energy level, and one weak transition from the H 13 2 CO isotopologue. The H 2 CO and H 13 2 CO are fitted simultaneously in XCLASS using an isotopic ratio calculated from Wilson & Rood (1994): 12 C/ 13 C ≈ 7.5 × d gal + 7.6. The 12 C/ 13 C ratio is listed in Table 1, for each region. While formaldehyde is a good low-temperature gas tracer at T kin < 100 K (Mangum & Wootten 1993), at high densities and temperatures, the 3 0,3 −2 0,2 transition is optically thick and a reliable temperature can no longer be derived from the line ratios with this method (Rodón et al. 2012;Gieser et al. 2019). To determine temperatures at higher densities, we use nine methyl cyanide lines (J = 12−11 and K = 0−8, E u /k B = 69−526 K). In addition, seven lines of the CH 13 3 CN isotopologue are fitted simultaneously (J = 12−11 and K = 0−6, E u /k B = 69−326 K). A detailed discussion of the line optical depth is given in Sect. 4.2.
In XCLASS, each molecule can be described by a number of emission and absorption components. The fit parameter set of each component consists of the source size θ source , the rotation temperature T rot , the column density N, the linewidth ∆v, and the velocity offset from the systemic velocity off . It is possible to choose a variety of algorithms that can also be combined in an algorithm chain in order to find the best-fit parameters by minimizing the χ 2 value. We adopted an algorithm chain with the Genetic and the Levenberg-Marquardt (LM) methods optimizing toward global and local minima, respectively.
For each region, we use the myXCLASSMapFit function to fit the H 2 CO and CH 3 CN lines in each pixel within the primary beam. All lines with a peak intensity >10σ line,map are fitted with one emission component. We chose a high threshold of 10σ line,map so multiple transitions have a high S/N to accurately determine the rotation temperature. Only a single value can be set as the threshold in the myXCLASSMapFit function, therefore, toward the edges of the primary beam, the temperature estimates are less reliable due to an increase in the noise.
Under the assumption of local thermal equilibrium (LTE), the kinetic temperature of the gas can be estimated from the rotation temperature T kin ≈ T rot . As HMSFRs have high densities n > 10 5 cm −3 toward the locations of the protostars, the LTE can be assumed here (Mangum & Shirley 2015). The critical densities, n crit = A ul C ul , for the CH 3 CN and H 2 CO lines are ∼4 × 10 6 cm −3 and ∼3 × 10 6 cm −3 , respectively (A ul is the Einstein A coefficient and C ul is the collisional rate coefficient). Here, we use C ul measured at 140 K taken from the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005). The H 2 CO and CH 3 CN temperature maps of each region are shown in Fig. C.1. The H 2 CO lines trace the extended low temperature gas at 10−100 K. The CH 3 CN maps mostly show spatially compact emission tracing the higher temperature gas at >100 K.
For each selected position (Table A.1), we extract the H 2 CO and CH 3 CN radial temperature profile which are binned in steps of half a beam (∆r). In each bin, the uncertainties are computed from the standard deviation of the mean. After a first visual inspection of all radial profiles of the 120 positions, we select all positions with a radially decreasing temperature profile along at least two beams in at least one temperature tracer. With this method, we are able to extract in total 22 positions, which we define as "cores". The radial temperature profiles of the cores A66, page 9 of 68 are shown in Fig. 4. In cases where both temperature tracers are detected toward the central core (e.g., IRAS 23033 2), the observed H 2 CO temperature profile has significantly higher values when compared to CH 3 CN. As discussed previously, this is due to the high optical depth in these dense regions causing the fit algorithm to converge toward high temperatures. The line optical depth is discussed further in Sect. 4.2. We fit the profiles with a power-law according to Eq. (1) using the minimum χ 2 method to derive the temperature powerlaw index q. We used the CH 3 CN temperature profile for the fit if it is detected along at least two beams and the H 2 CO temperature profile otherwise. The inner radius is the temperature at a radius of half the beam size r in = ∆r. The outer radius r out is chosen as a local minimum, when T rot (r out ) < T rot (r out + ∆r) and T rot (r out ) < T rot (r out + 2∆r). In cases where the outer radii of nearby cores overlapped or the 2D temperature distribution would become highly asymmetric, we reduced the outer radius. The outer radii are also marked in Fig. 4. The fit results for q, r in , r out , and T kin (r in ) are summarized in Table 3 for each core.
A histogram of the temperature power-law index q is shown in Fig. 5. Most of the data points are distributed between 0.1 and 0.6, but three outliers are located at q > 0.9. A Gaussian fit to the data with q = 0.1−0.6 yields an average value of q = 0.4 ± 0.1, which is in very good agreement with theoretical predictions (Emerson 1988;Osorio et al. 1999;van der Tak et al. 2000) and observations (e.g., Palau et al. 2014Palau et al. , 2020. The issue of optical depth of the H 2 CO lines and CH 3 CN being only detected toward the densest parts result in the high uncertainties and spread in q. We observe a broad range of q with shallow profiles (q = 0.1) up to steep profiles (q = 1.5). In order to study if q is constant for all high-mass cores or if the range of q is real, observations of more cores are required. Osorio et al. (1999) and van der Tak et al. (2000) suggest steeper values for q on scales <2000 au. In addition, the temperature maps and radial profiles are 2D projected, whereas the real 3D profile may be more complicated.

Density structure
In contrast to the merged (NOEMA + IRAM 30 m) spectral line data, the 1.37 mm continuum NOEMA data suffer from missing flux due to the lack of short-spacing information. Therefore, it is difficult to reliably derive radial intensity profiles in the final images that are required to determine the density structure. This issue can be minimized by analyzing instead the complex visibilities of the 1.37 mm continuum data in the uv-plane (Adams 1991;Looney et al. 2003;Beuther et al. 2007b;Zhang et al. 2009). Assuming spherical symmetry, the power-law index of the complex visibility profile α is related to the density powerlaw index of the radial density distribution p by (Looney et al. 2003): For each core position, the phase center is shifted to this location. The remaining cores within a region are fitted as a point source + circular Gaussian (to account both for the compact and the extended emission) and subtracted using the UV_FIT task in GILDAS. With this approach, the blending of nearby cores in the visibility profiles within a single region can be minimized. The azimuthally averaged complex visibilities, computed using the UVAMP task in MIRIAD (Sault et al. 1995), are binned in steps of 20 kλ. The radial visibility profiles are shown in Fig. 6. The uv distances range from ∼20−800 m corresponding to spatial scales of ∼10 3 −10 5 au at distances of a few kpc. We apply a powerlaw fit to the data in order to derive the power-law index α. The density power-law index p is calculated according to Eq. (3), taking into account the temperature power-law index q derived in Sect. 3.2. The results for α and p for all cores are summarized in Table 3.
Most of the visibility profiles ( Fig. 6) can be described by a single power-law. The higher scatter at large uv distances is due to the fact that less data is available at long baselines. For IRAS 23385, the visibility profiles of core 1 and 2 do not follow a simple power-law. Toward smaller spatial scales (<10 4 au), the profiles are steeper. This could be due to the fact that even though the contribution of nearby cores is minimized by subtracting a point source + circular Gaussian, in the case of IRAS 23385, it is not possible to clearly distinguish the contribution of both cores, which are embedded within a common dust envelope.
The histogram of the density index p is shown in Fig. 7. Using the observationally derived values of q, the power-law index ranges from 1.4−2.4. Fixing the temperature power-law index q to the mean value of 0.4 and calculating p with the results from the uv-analysis, the distribution of p gets narrower, peaking at p = 2.0−2.1. A Gaussian fit to this distribution yields a mean of p = 2.0 ± 0.2. The results of the derived physical structure of these 22 cores are discussed in detail in Sect. 6.1.

Molecular gas content
In this section, we analyze the chemical content of the molecular gas toward the 18 CORE HMSFRs by studying the molecular column densities derived toward the 120 positions. The column density of molecular hydrogen N(H 2 ) is derived from the 1.37 mm dust continuum emission (Sect. 4.1). In addition, molecular column densities are derived from the merged spectral line data with XCLASS by fitting the emission lines assuming  M core can be derived from the continuum intensity I ν assuming that the emission is optically thin (Hildebrand 1983): with a gas-to-dust mass ratio η = 150 (η =  Table 1.4 and 23.1 in Draine 2011, respectively), mean molecular weight µ = 2.8 (Kauffmann et al. 2008), hydrogen mass m H , beam solid angle Ω, dust opacity κ ν = 0.9 cm 2 g −1 for dust grains with a thin icy mantle at a gas density of 10 6 cm −3 at 1.3 mm (Ossenkopf & Henning 1994), the Planck function B ν (T ), and distance d. F ν is the integrated intensity of each core within the outer radius, r out . A66, page 13 of 68 A&A 648, A66 (2021) Notes. The density power-law index p is calculated according to Eq. (3). The core mass M core is calculated according to Eq. (5). The estimate of the chemical timescale τ chem is explained in detail in Sect. 5.2. a(b) = a × 10 b . Fit results for r in , r out , T kin (r in ), and q derived from either the H 2 CO or CH 3 CN radial temperature profiles (Fig. 4) are indicated by either a † or * , respectively. We use T kin for the temperature, T, in the Planck function taken from the thermometers H 2 CO and CH 3 CN, assuming LTE conditions (T kin ≈ T rot ). If the spectrum is extracted from a core position, T kin is taken from the radial temperature fit described in Sect. 3.2 with T = T kin (r in ) (see Table 3). If the spectrum is extracted from a position not corresponding to a core, the kinetic temperature is computed from T rot (CH 3 CN) if detected or from T rot (H 2 CO) otherwise. If there is no temperature tracer detected toward the position, we adopt a lower limit of T kin = 20 ± 10 K, as the lowest derived rotation temperatures range between 10−30 K.
In order to validate that the assumption of optically thin dust emission is valid, the continuum optical depth, τ cont ν , is computed for each position (for a derivation of the equation, see Appendix A in Frau et al. 2010): The kinetic temperature, T kin , molecular hydrogen column density, N(H 2 ), and continuum optical depth, τ cont ν , are listed in Table A.1 for all 120 positions. The uncertainties of N(H 2 ) and M core are calculated based on the assumption of Gaussian error propagation and include the uncertainty of the continuum intensity, with an estimated 20% flux calibration uncertainty ) and uncertainty of the derived rotation temperature, ∆T kin (listed in Table A.1). The optical depth, τ cont ν , is 1 toward most positions, so optically thin dust emission can be assumed here and the H 2 column density and core mass calculation provide reliable results. The only exceptions with τ cont ν > 1 are positions 1 and 2 of the W3 H2O region, which corresponds to the W3 OH UCHII region (see also Ahmadi et al. 2018). The molecular hydrogen column density N(H 2 ) has a mean value of 1.5 × 10 24 cm −2 , ranging from 2.7 × 10 22 cm −2 up to 2.8 × 10 25 cm −2 . The core masses, M core , are listed in Table 3. We find a mean core mass of 4.1 M in a range between 0.3 M and 11.6 M . As discussed previously, due to missing shortspacing information, both N(H 2 ) and M core should be considered as lower limits. However, the high core masses indicate that they are harboring protostars that will eventually end up as massive stars.

Spectral line modeling with XCLASS
We use the spectral line data of the CORE project to derive molecular column densities of 11 different species using the XCLASS software. A description of the XCLASS software is presented in Sect. 3.2. Using the myXCLASSFit function, individual spectra are fitted species by species with one emission component to derive the molecular column density, N.
A spectrum is extracted from the merged spectral line data for each position listed in Table A.1. The noise σ line in each spectrum is computed in a line-free range from 219.00 to 219.13 GHz. Compared to the average map noise σ line,map listed in Table 1, the noise in each spectrum σ line may have higher or lower noise values since the noise distribution is not uniform within the primary beam. The systemic velocity LSR is determined by fitting the C 18 O 2−1 transition and may differ from the global LSR listed in Table 1 due to velocity gradients in the region, hence, it allows us to employ a narrow parameter range for off and so fitting strong nearby emission lines is avoided. The systemic velocity LSR and noise σ line are listed in Table A.1 for each position.
All molecules and their corresponding transitions that are fitted with XCLASS are listed in Table B.1. Lines which are blended with transitions of other detected species (at a resolution of 3 km s −1 ) are also listed, but excluded from the fit. For most molecules, only the rotational ground-state level = 0 can be detected in our spectral setup. For SO 2 and HC 3 N vibrationally excited levels ( x > 0) are present and for CH 3 OH torsionally excited lines are detected ( t = 1). The following species for which we observe multiple isotopologues are fitted simultaneously: OCS and O 13 CS; SO 2 and 34 SO 2 ; H 2 CO, and H 13 2 CO; HC 3 N and HCC 13 CN; CH 3 CN and CH 13 3 CN. The isotopologue ratios are summarized in Table 1 and are calculated either from Wilson & Rood (1994) ( 12 C/ 13 C ≈ 7.5 × d gal + 7.6) or taken from references within ( 32 S/ 34 S ≈ 22). The uncertainties of these isotopic ratios are high due to a large scatter of the data points. However, we do not observe a sufficient number of strong isotopic lines to measure it more precisely. We did not fit 13 CO and C 18 O simultaneously because the 13 CO 2−1 line has a high optical depth (see Table B.1).
We use an algorithm chain with the Genetic and Levenberg-Marquardt (LM) methods and to estimate the uncertainties of the fit parameters, we used the MCMC error estimation algorithm afterward. The following criteria are applied to the fitted model spectrum and column density, N +∆N upp −∆N low , of each species in order to determine bad fits: (1) model spectra with a peak intensity < 3σ line ; (2) the upper error of the column density ∆N upp is converging toward high values (∆N upp > 10 × N); and 3) the lower error of the column density is not constrained ∆N low = 0 cm −2 . With these three criteria, weak and unresolved lines are automatically discarded and we use the best-fit value of N only as an upper limit.
The column densities, N, for all species fitted with XCLASS and their uncertainties derived with the MCMC error estimation algorithm are summarized in Tables E.1-E.3 for all positions. Histograms of the logarithmic column density distributions including N(H 2 ) are shown in Fig. 8, with the mean and standard variation of the column density noted in each panel (upper limits are not included). The logarithmic bin width is set to 0.5. Separate histograms of the core and non-core populations are shown in the same panels. However, as discussed in Sect. 2, for the spectral line data, we have short-spacing information, while we do not have this for the continuum data, hence we systematically underestimate the H 2 column density.
There are species with a distribution having a clear column density peak (e.g., H 2 , 13 CO, C 18 O, SO, DCN, H 2 CO, HNCO). Other species have a double-peaked distribution with a clear separation between the core and the remaining positions (e.g., OCS, HC 3 N, CH 3 OH, CH 3 OH; t = 1, and CH 3 CN). In these cases, the column density is enhanced by a factor of ∼10−100 toward the core positions. There are not enough data points for the SO 2 , HC 3 N; 7 = 1, HC 3 N; 7 = 2 to draw any conclusions about the distribution, however, they are detected mostly in the densest regions toward core positions. In general, high column densities are found toward core positions and low column densities are found toward the remaining positions.
To account for the fact that toward the core positions the column density is expected to be higher in a higher density region, we show relative abundance N(X)/N(C 18 O) histograms in Fig. D.1 (upper limits are not included). Assuming that both species trace the same emission region, the relative abundances are independent of the absolute column density value, which differ from region to region and from core to core. Normalized to N(C 18 O), most species have a single-peaked distribution. Exceptions are OCS, HNCO, HC 3 N, CH 3 OH, and CH 3 CN with double-peaked distributions indicating that high temperature gas-phase chemistry has a big impact on these N-bearing molecules by increasing their abundances. However, in general there is still a clear difference between core positions (high abundance) and non-core positions (low abundance). The fact that larger molecules have a clearer distinction between the core and non-core positions, while for simpler species it is less obvious (e.g., 13 CO and DCN), hints that the emission of COMs is associated with the cores while simple molecules are abundant in envelope as well. The difference of the core and the remaining positions indicates that the high densities and possible energetic processes around the protostars have a strong impact on the molecular abundances in the gas phase (e.g. through protostellar outflows, shocks, disk accretion, and strong radiation from the protostars). Correlations of the derived column densities are discussed in Sect. 6.3.
Observed and XCLASS modeled spectra are shown in  observed spectra vary from region to region, but also within a region. Typical hot core spectra are observed for the positions AFGL 2591 1, CepA HW2 1, G75.78 1, W3 H2O 1, and W3 H2O 2. Many weak emission features are detected in the spectra at a ∼2−3σ line level originating from COMs. These COM emission features are difficult to fit as the transitions are weak, have similar upper energy levels E u /k B , and are blended at a spectral resolution of 3 km s −1 , so they were excluded from this analysis. A detailed study of the integrated line emission (including line stacking of weak COM emission lines) will be presented in a future study (Gieser et al. in prep.). Fewer species are detected in spectra toward the non-core positions. In contrast to the line-rich sources, there are several sources that show only a handful of emission lines mainly from CO-isotopologues, SO, and H 2 CO even toward the continuum peak positions. The linepoor regions are G139.9091, G138.2957, S87 IRS1, and S106. A66, page 18 of 68 In the case of IRAS 23033, core 1 has a line-poor spectrum, whereas for core 2 and 3, which are embedded in a common envelope, significantly more emission lines are detected. These cores either could be at different evolutionary stages or they are embedded in an inhomogeneous local radiation field. The former can be investigated by applying a chemical model to the observed column densities and by estimating the chemical ages of the regions. We investigate this possibility in the next section.
The mean and maximum line optical depth, τ line mean and τ line max , computed for each fitted transition with XCLASS are summarized in Table B.1. The 13 CO 2−1 transition has the highest mean optical depth with τ line mean = 1.5, but for most other species and transitions, the mean line optical depth is <1, so the column density and temperature determination should be reliable, especially when many transitions of a molecule are fitted simultaneously. The mean optical depth of the H 2 CO 3 0,3 −2 0,2 transition is a factor of four higher than the remaining two transitions. When the 3 0,3 −2 0,2 transition is optically thick, temperature estimates depending on the line ratios are difficult to determine with the remaining two optically thin lines as they have similar upper energy levels (68 K). In 23 out of the 118 spectra where H 2 CO was detected and fitted with XCLASS, the calculated optical depth of the 3 0,3 −2 0,2 transition is τ line > 1.

Physical-chemical modeling of the cores
The continuum data of the CORE sample show a large diversity in fragmentation properties  and our analysis (shown in Sect. 4) demonstrates that the composition of the molecular gas varies within as well as between the regions: some have a rich plethora of molecular lines, while others have line-poor spectra. This diversity of physical and chemical properties could be explained by a number of reasons. Magnetic fields or different initial density structures could explain the variety in fragmentation properties . Different initial conditions in, for instance, the large-scale kinematics and mass distribution might also have an effect on the molecular abundances. To investigate whether the observed variation of the physical and chemical properties of the cores may be due to the fact that the CORE regions have different ages, we model the chemical evolution of the 22 cores in the following.

MUSCLE setup
We applied the physical-chemical model to the physical properties and molecular column densities of each core determined from the CORE 1 mm observations in order to estimate the chemical ages. MUSCLE (MUlti Stage ChemicaL codE) has previously been successfully applied to the CORE pilot regions NGC 7538 S and NGC 7538 IRS1 (Feng et al. 2016) as well as the CORE region AFGL 2591 (Gieser et al. 2019). The model comprises spherically symmetric physical structures. The temperature and density profiles of the cores are described by power-laws up to the outer radius, r out , with index q and p, respectively (see Eqs. (1) and (2)). At an inner radius, r in , and further in, the density and temperature reach a constant value. We adopt 40 logarithmic grid points for the radial profiles.
On top of this static physical structure, the time-dependent gas-grain chemical network ALCHEMIC (Semenov et al. 2010) computes the abundances of hundreds of atomic and molecular species using thousands of reactions. A detailed description of MUSCLE can be found in Gerner et al. (2014Gerner et al. ( , 2015. We adopt most of the model parameters from the AFGL 2591 case study described in Gieser et al. (2019), which yielded a good estimate Initial Chemical Abundances HMPO model Best-fit IRDC stage (g) Inner radius r in 12 700 au Outer radius r out 0.5 pc Temperature T in at r in 11.3 K Temperature power-law index q 0.0 (isothermal) Density n in at r in 1.4(5) cm −3 Density power-law index p 1.5 Stage lifetime τ IRDC 16 500 yr HMC model Best-fit HMPO stage (g) Inner radius r in 103 au Outer radius r out 0.5 pc Temperature T in at r in 75.8 K Temperature power-law index q 0.4 Density n in at r in 1.5(9) cm −3 Density power-law index p 1.8 Stage lifetime τ HMPO 32 000 yr UCHII model Best-fit HMC stage (g) Inner radius r in 1 140 au Outer radius r out 0.5 pc Temperature T in at r in 162.9 K Temperature power-law index q 0.4 Density n in at r in 1.3 (8)  of the chemical age of this hot core compared to literature estimates. A summary of the input parameters is listed in Table 4. In contrast to the AFGL 2591 case-study, we use a higher value for the cosmic ionization rate ζ CR based on a study of multiple HMSFRs by Indriolo et al. (2015). These authors find that ζ CR is constant at a Galactic radius > 5 kpc and with all the CORE regions at Galactic distances >7 kpc (Table 1), we use a constant value of 1.8 × 10 −16 s −1 . By setting the extinction at r out to A v = 10 mag , the core is shielded from the interstellar ultraviolet radiation field. For each of the 22 cores, we run a physical-chemical model with MUSCLE. The input are the H 2 column density (Sect. 4.1) and all molecular column densities derived with XCLASS (Sect. 4.2). The CO column density N(CO) is calculated from N(C 18 O) as C 18 O is less optically thick than 13 CO (  Table 1 for each region. For HC 3 N and CH 3 OH, we compute the mean column density of the rotational ground state and vibrationally or torsionally excited states for the MUSCLE input. Molecular column densities, for which only upper limits could be determined, are also set as upper limits in MUSCLE. We set the temperature structure of the model core to the observed temperature profile (Sect. 3.2) and used the density power-law index derived from the continuum visibility analysis (Sect. 3.3).
Two undetermined model parameters remain. We do not know, for one, how evolved the cores are, described by the chemical age, τ chem , and second, what the initial chemical composition of the parental molecular cloud/clump was. Due to the fact that the CORE regions are far more evolved than typical cold IRDCs and the physical structure of each model stage is static, we have to define a sensible initial condition for the chemical composition. The initial conditions we apply are based on a study of 59 HMSFRs using single-dish observations (Gerner et al. 2014(Gerner et al. , 2015. These HMSFRs were classified according to their evolutionary stage (IRDCs, HMPOs, HMCs, and UCHII regions) and a template was created from the average column densities for each evolutionary stage. The four template stages were modeled using MUSCLE to create average abundances for all molecular and atomic species and to estimate a mean chemical age of each evolutionary stage. The properties of their template IRDC, HMPO, and HMC model are summarized in Table 4. Following the convention by Gerner et al. (2014Gerner et al. ( , 2015, the chemical age τ chem is 0 yr when the gas density reaches 10 4 cm −3 . Based on the temperature profiles of the 22 cores, we can assume that they lie somewhere between the HMPO and early UCHII stage, since the average temperatures around cores are too high to be classified as IRDCs (T ∼ 20 K). There are a few known UCHII regions with strong free-free emission at cm wavelengths resolved in the CORE data. In W3 H2O, the Western part (around position 1 and 2 in Fig. 3) is the UCHII region W3 OH. In W3 IRS4, the Southern ring-like structure is an UCHII region as well (Mottram et al. 2020). However, for the W3 OH UCHII region we do not find a clear radial decreasing temperature profile and toward the W3 IRS4 UCHII region no H 2 CO or CH 3 CN line emission is detected at a 10σ line,map level to estimate the kinetic temperature. The S106 region is a UCHII as well, for which we do not detect neither H 2 CO nor CH 3 CN emission around the compact core. This already suggests that toward this later stage, the molecular richness in these regions is decreased.
To test which initial conditions (see Table 4) fit best to the observed molecular column densities, we model each of the 22 cores with initial abundances after an initial IRDC phase (referred to as the HMPO model), after an initial HMPO phase (referred to as the HMC model), and after an initial HMC phase (referred to as the UCHII model). While most cores are unlikely to have formed a strong UCHII region yet, we include the UCHII model, as the observations in Gerner et al. (2014Gerner et al. ( , 2015 have large beam sizes (11 and 29 ), and UCHII regions may have contamination from less evolved line-rich objects, which are blended into the single-pointing spectra. It is not our aim to classify the cores into these evolutionary stages, but to find sensible initial chemical conditions as an input for our physical-chemical modeling. For example, an evolved HMC that is more evolved than the template HMC from Gerner et al. (2015) would best be described best by the UCHII model in our nomenclature.
For each model, the chemical evolution, τ model , runs up to 100 000 yr in 100 logarithmic time steps. In each time step, the computed radial abundance profiles are converted into beamconvolved column densities with the beam size fixed to the mean synthesized beam of the observations. The best-fit model is determined by a minimum χ 2 analysis by comparing the modeled and observed column densities in each time step and for all three adopted initial condition models. Applying this physical-chemical model allows us to estimate the chemical age.

Chemical ages
The best-fit chemical age τ chem , χ 2 , and percentage of wellmodeled molecules are shown in Table G.1 for each initial abundance model and core. Gerner et al. (2014) estimated that chemical ages are uncertain by a factor of between two and three and that the modeled column densities are uncertain by a factor of ten. But this depends on the number of modeled molecules, but also cores embedded in complex dynamic environments are harder to fit with our model. The chemical age τ chem is the sum of the time of the initial abundance model and τ model , so for the HMPO model: τ chem = τ IRDC + τ model , for the HMC model: τ chem = τ IRDC + τ HMPO + τ model ; and for the UCHII model: τ chem = τ IRDC + τ HMPO + τ HMC + τ model .
For some cores, multiple initial condition models have a similarly low χ 2 (e.g., the HMPO and UCHII model for core 1 in IRAS 23033). In these cases, we cannot constrain the chemical age well. Comparing the lowest χ 2 model with the remaining initial condition models, if the χ 2 difference is less than 5%, only chemical age ranges spanning over these models are further considered. Table 3 shows the chemical age τ chem for models with a clear lowest χ 2 initial condition model or a time range in chemical age for cores with multiple best-fit initial condition models. The estimated chemical timescales τ chem of the cores vary between ∼20 000 and 100 000 yr within the regions of the CORE sample with a mean of ∼60 000 yr. The youngest core being G084.9505 1 and the oldest core being G108.75 2.
A comparison of the best-fit modeled and observed column densities for all cores is shown in Fig. 9. For most cores, the model underestimates the H 2 CO and CH 3 OH column densities as compared to the observed values. This can be explained by the fact that the quasi-static model does not sufficiently take into account the warm-up stage from 30−80 K where surface chemistry on the dust grains plays an important role and where these two molecules are formed by subsequent hydrogenation of CO. These discrepancies between modeled and observed H 2 CO and CH 3 OH column densities have already been noticed by Gerner et al. (2014) in their template HMPO stage modeled with MUSCLE. They explain that this is due to the fact that the formation route of H 2 CO consists of grain-surface as well as gas-phase chemistry which are both time-dependent and not correctly implemented in the chemical models. This results in the over and underproduction of these species, which is also the case in our modeling results, shown in Fig. 9.
Large discrepancies between the modeled and observed column density exist also for the SO molecule for which the model overproduces the SO column density by a factor >10 for many cores (e.g., IRAS 23033 core 1, 2, and 3). This overproduction in SO is seen in all three initial condition models, but in most cases, other modeled S-bearing species (OCS and SO 2 ) are modeled well. The applied initial chemical conditions based on the Gerner et al. (2015) models also included S-bearing species (SO, CS, and OCS). Their initial IRDC stage model started with elemental abundances and only H 2 in molecular form taken from the low metals set of Lee et al. (1998). But in order to fit the IRDC phase accurately, Gerner et al. (2015) had to increase the initial A66, page 20 of 68  However, an overproduction of SO is also seen in their best-fit HMPO, HMC and UCHII models. This might be connected to a poorly understood chemistry of the reactive SO molecule, as also in their models the remaining S-bearing species can be reproduced properly. In addition, as only one SO transition is covered in our spectral setup, which can typically be optically thick (Table B.1 and Fig. F.1), we may underestimate the observed SO column density. This might partially explain the differences between the modeled and observed SO column density. With multiple cores resolved within a region, it is possible to study how the chemical timescale τ chem varies across small spatial scales. In the IRAS 23033 region, core 1 seems to be more evolved (∼30 000−100 000 yr), even though the spectrum is line-poor (see Fig. F.1) compared to the spectra of core 2 and 3 which are embedded in a common envelope (see A66 Fig. 3) and for which we estimate similar chemical timescales of ∼20 000 yr. Core 1 and 2 in CepA HW2 have a chemical age of ∼80 000 yr and ∼20 000−90 000 yr, respectively. The cores are very close (∼2300 au), but within our sensitivity limit, these cores are not embedded in a common envelope, but have very steep density profiles (p 2). In IRAS 21078, core 1 and 2 have a chemical age of ∼60 000−90 000 yr and ∼50 000 yr, respectively, suggesting a small age gradient. The cores are embedded within a common envelope and have small projected separations. Core 3 and 4 in W3 H2O have a chemical age of ∼90 000 yr and ∼20 000−90 000 yr, respectively. In the IRAS 23385 region, core 1 is younger (∼50 000 yr), while core 2 is estimated to be much older (∼100 000 yr). In G108.75 a large difference between the chemical ages of core 1 (∼20 000 yr) and core 2 (∼110 000 yr) is estimated. The cores have a separation of ∼20 000 au, but have the same systemic velocity (Table A. radiation field or complex dynamics could be the reason for this large chemical age difference. One of the limitations of MUSCLE is that the physical structure (radial temperature and density profiles) is static within each evolutionary stage (IRDC, HMPO, HMC, and UCHII). In reality, these properties do change on timescales smaller than the chemical timescales derived here; in addition, the dynamics (e.g., gas inflow) are important factors to consider. Currently, 3D timedependent physical models in combination with a full chemical network are computationally expensive. Therefore, we use the approach of our quasi-static physical model by considering the four different evolutionary stages. More sophisticated physicalchemical models in the future are required to include 3D gas dynamics and the evolution of the density and temperature structure. In addition, a larger number of molecular column densities would better constrain the model parameter space.

Physical structure of high-mass star-forming cores
Various methods were applied in the literature to observationally derive the density profiles of envelopes in HMSFRs (e.g., summarized in Table 6 in Gieser et al. 2019). Some of these studies are based on observations with single-dish telescopes with beam sizes >10 tracing the clump-scale envelope, while interferometric observations trace the core-scale envelope. We selected studies in the literature for which the density structure was determined for a sample of cores or clumps and we extracted the typical sizes r from their studies. The results in comparison with our study are shown in Fig. 10. At scales of 1 pc down to 0.01 pc, it seems that the density index p lies between 1.7 and 2.0, which is close to the values inferred in low-mass star-forming regions (Motte & André 2001). To investigate this further, we observed the CORE regions with the NIKA2 instrument at the IRAM 30 m telescope and an analysis of the density structure at clump-scales will follow (Beuther et al., in prep.).
The observationally-derived density and temperature profiles (q = 0.4 ± 0.1 and p = 2.0 ± 0.2) are in agreement with theoretical studies of HMSF, but the physics of how massive stars form is still not fully understood. Currently, theoretical models propose the formation of high-mass stars through: a monolithic collapse of turbulent cores (McKee & Tan 2002; protostellar collisions and coalescence in dense clusters (Bonnell et al. 1998;Bonnell & Bate 2002); or competitive accretion in clusters (Bonnell et al. 2001;Smith et al. 2009;Hartmann et al. 2012;Murray & Chang 2012). The density and temperature structure are important parameters of the initial cloud and proceeding clump and core scales. For example, early star formation models by Shu (1977) and Shu et al. (1987) that model the gravitational collapse of an isothermal sphere find that p = 2 in the outer envelope and p = 1.5 in the inner region where the gas is free-falling onto the central region. McLaughlin & Pudritz (1996 used a logatropic equation of state and a non-isothermal sphere and find that at an initial density profile of p = 1, the profile steepens to p = 1.5 after the collapse. In the turbulent core model by McKee & Tan (2002 the authors assume p = 1.5 based on observational constraints. In Bonnell et al. (1998) the density profile in the outer region has the form p = 2 and a shallower, near-uniform profile in the central region. Murray & Chang (2012) explore their models by varying p from 0 (uniform), 1, and 2 (isothermal).
The density structure is important for the physical and chemical evolution of HMSFRs. It is therefore important to quantify the initial density profile on cloud to clump and cores scales and how it changes with time. While theoretical models usually do not predict, but rather assume a given density profile, observations of HMSFRs on different scales can help to narrow down the parameter space (see Fig. 10). Hydrodynamic simulations reported by Chen et al. (2021) investigate how changes of p in giant molecular clouds with an initial radius of 20 pc affect massive star cluster formation. They find that for steep density profiles, p = 2, there is a centrally concentrated cluster, while for shallower profiles, hierarchical fragmentation occurs. Hydrodynamic simulations from Girichidis et al. (2011) show that massive protostars form only in clouds with a density index of p = 1.5 or p = 2.0, while for uniform or Bonnor-Ebert-like profiles, a large fraction of low-mass stars form. The authors found that turbulence and the initial density profile are important aspects for the evolution of the cloud and the formation of clusters.
We study the correlations of all core properties shown in Table 3 using the Spearman correlation coefficient r. This statistical tool can be used to check if two data sets have a positive correlation (r = 1), negative correlation (r = −1), or no correlation (r = 0). We assume that a high correlation exists if r > 0.8. For example, Feng et al. (2020) finds a negative correlation for the H 2 column density and dust temperature for cold highmass clumps using the Spearman correlation coefficient, r. A big advantage compared to the Pearson correlation coefficient is that linear, as well as nonlinear, correlations are considered in the calculation of r. In addition, we add the M/L ratio of the region listed in Table 1 in Beuther et al. (2018) as a parameter for each core. However, the interpretation is difficult as multiple cores within a region have the same M/L ratio. A mean chemical age is used in the computation of r for cores for which only a time range can be estimated (Table 3).
The results for the correlation coefficient, r, are shown in Fig. 11, where all pairs with a correlation >0.8 are highlighted.
Unfortunately, a small sample of only 22 cores does not allow us to find many strong correlations. Observations of many HMS-FRs at core-scales are required to study these relations in a better statistical way, which will be possible, for example, with the ALMAGAL survey, an ongoing ALMA large program observing  Fig. 11. Spearman correlation coefficient r for the derived physical and chemical core parameters listed in Table 3. Values higher than 0.8 are marked by a green circle. more than 1 000 HMSFRs. A high correlation is found between the inner and outer radius, which is due to the fact that we are resolution-limited and the regions are located at different distances. The correlation between p and α is due to Eq. (3). However, we find a strong negative correlation between the M/L ratio and the H 2 column density of the cores, so more evolved cores have a higher beam-convolved H 2 column density. The M/L ratio, proposed to be a good tracer of the evolutionary stage, is investigated in the following section. Sridharan et al. (2002) found that the distance-independent ratio M/L of UCHII regions is lower (∼0.005) than the ratio of HMPOs (∼0.025), as UCHII regions are more evolved and thus more luminous. This has also been confirmed by observations of large samples of HMSFRs, for instance, by Molinari et al. (2008Molinari et al. ( , 2010Molinari et al. ( , 2016Molinari et al. ( , 2019; Maud et al. (2015); Urquhart et al. (2018). In Fig. 12 we plot the region-average M/L ratio taken from Table 1 in Beuther et al. (2018) against the estimated chemical ages τ chem . The chemical timescale and a factor 2 uncertainty is shown for cores with a clear best-fit model. For cores with estimated chemical age ranges, the mean chemical age are shown with error bars spanning over the time range. The luminosities L have uncertainties on the order of ∼30% (Mottram et al. 2011b) and the masses M calculated in Beuther et al. (2018) are expected to be uncertain within a factor of 3. We also show the average M/L ratio for HMPOs and UCHII regions derived by Sridharan et al. (2002). We corrected the M/L ratios of Sridharan et al. (2002) by a factor of 0.5, as M was taken from Beuther et al. (2002) where the reported values for M were lower by a factor of 2 (Beuther et al. 2005). The cores G108.75 1 and 2 are excluded here, because no consistent continuum data is available to reliably derive the mass ). There is a tendency such that the older cores have a lower M/L ratio. However, as the M/L ratios are determined on much larger scales and we resolve multiple cores within each region, it is difficult to compare these properties.

M/L ratio as a tracer of evolutionary trends
Estimating the chemical ages with MUSCLE shows that a line-rich spectrum does not have to imply that a core is more evolved compared to a core with a line-poor spectrum. For example, the line-poor core 1 in IRAS 23033 is estimated to be older than the line-richer cores 2 and 3. When a hot core evolves to become an UCHII region, the strong protostellar radiation destroys molecules. Hence line-poor spectra can be found at early and at late evolutionary stages. Observations of more species would be helpful to better constrain the chemical model.
We are not able to derive radial temperature profiles from the temperature maps of G138.2957, G139.9091, and S106 ( Fig. C.1). In G139.9091 and S106 there is no emission at the adopted 10σ line,map level. There is diffuse emission of H 2 CO in G138.2957, but it is too diffuse to derive a radial temperature profile. The spectra of these regions are line-poor and only simple species are detected ( 13 CO, C 18 O, SO, H 2 CO, see Fig. F.1). Based on the shape of the 1.37 mm continuum emission, G139.9091 and S106 have isolated cores with no significant envelope emission, while G138.2957 has diffuse dust emission (Fig. 3). Therefore one may conclude that G139.9091 and S106 are already more evolved and probably have τ chem > 100 000 yr. S106 is a well studied bipolar HII region (e.g., Roberts et al. 1995;Schneider et al. 2007). G139.9091 is associated with a HII region as well (e.g., Kurtz et al. 1994;Purser 2017;Obonyo et al. 2019). G138.2957 could be in a very young strongly depleted phase with τ chem < 20 000 yr or could be an evolved region and the diffuse emission is due to the disruption by the protostellar radiation. With observations of G138.2957 at 5.8 cm and 20 cm, a core component with an associated synchrotron jet is seen toward the location of the 1.37 mm continuum peak and the position has an associated infrared source (Obonyo et al. 2019). This suggests that G138.2957 is an evolved embedded cluster and not a young region. Deep observations at radio wavelengths (∼5−50 GHz) would also provide information on the presence of UCHII regions within the CORE sample. We find that for the known HII regions, we only detect simple species such as CO isotopologues, H 2 CO suggesting that a large fraction of the larger molecules are already destroyed in this stage.

Correlations between chemical species
In this section, we aim to establish which molecules show a correlation via chemical links or temperature effects. We compute the Spearman correlation coefficient r for the molecular column density pair combinations relative to N(C 18 O). SO 2 , HC 3 N; 7 = 1, and HC 3 N; 7 = 2 are excluded from this analysis as there are fewer than 20 column density data points. Ideally, we would compare the correlations of relative abundance pairs relative to N(H 2 ), but due to the issue of the missing flux, we use N(C 18 O), which is also detected toward all positions and is optically thin (Table B.1). Comparing the correlation between column density pairs, we ultimately face the issue that toward higher densities, the column density is also higher (as discussed in Sect. 4.2). The core and non-core positions were equally considered in the calculation of r and the results for all column density pairs are shown in Fig. 13.
In general, all pairs show a positive correlation, r > 0. High correlation coefficients (r > 0.8) are found between pairs of the following molecules: HC 3 N-SO, HC 3 N-OCS, CH 3 CN-SO, CH 3 CN-OCS, CH 3 CN-H 2 CO, CH 3 CN-HNCO, CH 3 CN-HC 3 N, CH 3 OH-CH 3 OH; t = 1, CH 3 CN-CH 3 OH. The lowest correlations are found for 13 CO and DCN, where there is no correlation with any species. In the case of 13 CO this is most likely due to a high optical depth (Table B.1), so the column density cannot be reliably derived. The case of DCN is more puzzling, but a more detailed study of the deuteration would be required. Unfortunately, the CORE spectral setup covers only the DCN 3−2 line, so follow-up observations of deuterated species are necessary to study how deuterium chemistry behaves on such scales.
Methyl cyanide (CH 3 CN) shows a strong positive correlation with most other species, including S-bearing and N-bearing species. A correlation of HC 3 N-CH 3 CN can be explained by gasphase N-chemistry in the envelope gas (Bergner et al. 2017). We find that OCS is correlated with other dense gas tracers (HC 3 N and CH 3 CN). A correlation exists between CH 3 OH and CH 3 CN even though there is no chemical link between methanol and methyl cyanide. Such a correlation has also been found toward low-mass star-forming regions (Bergner et al. 2017;Belloche et al. 2020) and toward the massive star-forming region G10.6−0.4 (Law et al. 2021). Belloche et al. (2020) argued that this is a temperature effect caused by chemically unrelated species having been evaporated from icy grain mantles by energetic processes. Additional high angular resolution observations of SiO would be helpful for studying the impact of shocks in greater detail. In the 1D physical-chemical modeling of HMS-FRs with MUSCLE by Gerner et al. (2014), CH 3 CN, and CH 3 OH are co-spatial in radial abundance profiles toward the inner hot core region. Urquhart et al. (2019) find that the line integrated ratios of S-and N-bearing species are positively correlated with the dust temperature in a sample of high-mass star-forming clumps. Gratier et al. (2013) found that CH 3 CN is much more abundant in the photo-dissociation region (PDR) of the Horsehead nebula than the associated cold and dense core. This is consistent with the detection by Purcell et al. (2006) of 3 mm lines of CH 3 CN in 58 candidate hot molecular cores in a sample of 83 methanol maser-selected star-forming regions. The authors detected CH 3 CN in isolated methanol maser sites and found that CH 3 CN is more prevalent and brighter when an UCHII region is present, independently of the distance to the source. Guzmán et al. (2014) proposed that correlated abundances of CH 3 OH and CH 3 CN could be related to photochemistry, for instance, by photodesorption.
The column densities relative to N(C 18 O) show in general positive correlations. The N-bearing and S-bearing species seem to be chemically related by high temperature gas-phase chemistry. The strong non-correlation of DCN with any observed species requires high angular resolution follow-up observations of deuterated molecules toward these HMSFRs (e.g., N 2 D + , DCO + ). The correlation between CH 3 CN and CH 3 OH can be due to a mutual evaporation temperature or due to photodesorption.
Computing the Spearman correlation coefficient r of the chemical timescales, τ chem , (Table 3) with the observed column densities N (Fig. 8) and with the relative abundances N/N(C 18 O) (Fig. D.1) for all 22 cores, we do not find strong correlations for most of the molecular column densities and relative abundances with the chemical age. We use the mean chemical age for cores for which only a time range can be estimated. Positive correlations of the column density and relative abundance with the chemical timescale (r = 0.5−0.9) are found for SO 2 , HNCO, HC 3 N; 7 = 1, and CH 3 OH; t = 1. We only detect HC 3 N; 7 = 2 in cores, which have a chemical age estimated to be >80 000 yr. This suggests that vibrationally and torsionally excited states of molecules are good indicators for a more evolved region, as these transitions have higher upper energy levels (Table B.1) and are only excited at a high kinetic temperature. We propose that SO 2 and HNCO are also good tracers of the evolutionary stage, however, this has to be investigated using a larger statistical sample.

Comparison of physical and chemical timescales
In the following, we compare the estimated chemical timescales with commonly applied physical timescales. The free-fall timescale, τ ff , is the time it takes for a spherical object to fully collapse under the influence of gravity with no additional forces. It only depends on the initial density, ρ: where G is the gravitational constant and ρ(H 2 ) the initial H 2 mass density.
The crossing timescale is the time it takes to cross the system once: with clump radius, R, and velocity dispersion, .
A comparison between the free-fall timescales for corresponding H 2 number densities of 10 5 cm −3 , 10 6 cm −3 , 10 7 cm −3 , and 10 8 cm −3 , crossing timescale (assuming R = 1 pc and = 4 km s −1 , Elmegreen 2000), and the chemical ages of the regions estimated with MUSCLE is shown in Fig. 14. A mean chemical age is used for cores for which only a time range can be determined. We then calculate the mean chemical age for regions associated with multiple cores. In these cases, the errorbars indicate the range of chemical ages. For regions associated with only one core, we assume an uncertainty of a factor of two for cores with a clear best-fit chemical age or the errorbars indicate the estimated time range for cores for which we can only determine a time range. The derived chemical ages are in agreement with an initial clump density of 10 5 −10 6 cm −3 . Indeed, observations suggest that IRDCs have densities n(H 2 ) > 10 5 cm −3 (e.g., Carey et al. 1998). The estimated crossing time is an upper limit, in agreement with the scenario that suggests that star formation occurs within a few crossing times (Elmegreen 2000). The agreement of the physical timescales with the chemical ages shows that even though the physical-chemical model is pseudotime-dependent, the estimates are feasible, especially when a homogeneous data set and analysis is used.
Other forces such as magnetic fields and turbulence can slow down or prevent collapse and are important in the inital diffuse ISM (Ballesteros-Paredes et al. 2007;Burge et al. 2016). However, in the diffuse ISM, molecules, with the exception of H 2 , have not formed yet. Our zero point for the chemical timescales is when the density reaches 10 4 cm −3 and for the freefall timescale, we also assume densities 10 4 cm −3 . In Beuther et al. (2018), we showed that the fragmentation scales that we observe are in agreement with thermal Jeans fragmentation, thus turbulent Jeans fragmentation is not significant here. This is fully consistent with the results found in a different sample by Palau et al. (2015). The role of magnetic fields in the CORE regions are being investigated by polarization observations with the Submillimeter Array (SMA).

Conclusions
In this work, we study the physical and chemical structure of 18 high-mass star-forming regions using the CORE 1 mm continuum and spectral line data. We quantify the chemical content for in total 120 positions including in total 22 cores. Combining the CORE observations with the physical-chemical model MUSCLE, we estimate the chemical timescales of the cores. Our main conclusions are summarized as follows: 1. Using H 2 CO and CH 3 CN line emission, we derived temperature maps of the regions. We identify "cores" as objects having a clear radially decreasing temperature profile. We fit these radial temperature profiles and obtained an average power-law index q = 0.4 ± 0.1, excluding three outliers with q ≥ 1, which is in agreement with theoretical predictions and calculations. 2. The 1.37 mm continuum visibility profiles of the cores reveal a mean density power-law index of p = 2.0 ± 0.2. Comparing these high-resolution density profiles to previous singledish and interferometric studies, the density profiles appear to stay roughly constant between p = 1.7−2.0 from scales of 1 pc down to 1000 au. The molecular hydrogen column density N(H 2 ) has a mean of 1.5 × 10 24 cm −2 toward the 120 positions and the core mass M core of the 22 cores has a mean of 4.1 M . 3. We derived molecular column densities of 11 species by fitting their spectral lines with XCLASS. Spearman correlation coefficients are evaluated for all molecule pairs. We find high correlations between N-bearing and S-bearing species that are chemically related through high temperature gas-phase chemistry. High correlations are also found for molecules that are not chemically related (CH 3 CN and CH 3 OH), but for which the correlation results from a common evaporation temperature or photo-desorption. 4. We applied the physical-chemical model MUSCLE to the observed column densities of each core to estimate the chemical age τ chem . We find a spread in age from τ chem ∼ 20 000−100 000 yr and a mean chemical age of ∼60 000 yr. Multiple cores within a region show that there can be an age gradient for largely separated cores (e.g., toward IRAS 23033 and G108.75), while close cores have a similar chemical age (e.g., toward IRAS 21078 and W3 H2O) suggesting a sequential star formation. 5. A strong correlation is found between the peak column density N(H 2 ) of the 22 cores and the M/L ratio of the region. In addition, we find a trend which posits that older cores have a lower M/L ratio. However, a larger sample is required to study the physical and chemical properties on core-scales in more detail. 6. We compare the chemical age with physical timescales.
We find that the chemical age is consistent with a free-fall timescale of the initial clump, at a density of 10 5 cm −3 and 10 6 cm −3 , which is consistent with density values typically found toward IRDCs. The chemical ages are smaller than a crossing time of the parental clump.
A66, page 26 of 68 7. We improved the quality (both intensity and noise) of the NOEMA continuum and spectral line data by applying selfcalibration. All of the standard and self-calibrated CORE data are publically available 4 . The CORE data reveal a large physical and chemical diversity on scales down to ∼300 au. Here, we use the molecular column densities of 120 positions toward the 18 high-mass star-forming regions to quantify the chemical content, with an emphasis on 22 cores. The case study of AFGL 2591 had already revealed complex structures around a single hot core (Gieser et al. 2019) using high-resolution interferometric observations. In a next step, we will study the spatial morphology of the molecular emission in order to evaluate the spatial correlations of the species. Notes. The entries are taken from the CDMS and JPL catalogs (Müller et al. 2005;Pickett et al. 1998).  Table 2 for values of σ cont for each region). Each core is marked in red and the dashed red circle indicates the outer radius of the radial temperature fit (Sect. 3.2). The beam size is shown in the bottom left corner in each panel. The pink bar in the top left corner indicates a linear spatial scale of 10 000 au. The primary beam size is indicated by a black circle and for regions with no extended H 2 CO temperature map a smaller field of view is shown.