# Kinematics and stability of high-mass protostellar disk candidates at sub-arcsecond resolution

A. Ahmadi, H. Beuther, F. Bosco, C. Gieser, S. Suri, J. Mottram, R. Kuiper, T. Henning, Á. Sánchez-Monge, H. Linz, et al.

► To cite this version:

A. Ahmadi, H. Beuther, F. Bosco, C. Gieser, S. Suri, et al. Kinematics and stability of high-mass protostellar disk candidates at sub-arcsecond resolution. Astronomy & Astrophysics - A&A, 2023, 677, pp.A171. 10.1051/0004-6361/202245580. hal-04295049

HAL Id: hal-04295049

https://hal.science/hal-04295049v1

Submitted on 2 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Kinematics and stability of high-mass protostellar disk candidates at sub-arcsecond resolution

## Insights from the IRAM NOEMA large programme CORE ★,★★

A. Ahmadi<sup>1,2 Ⓓ</sup>, H. Beuther<sup>2 Ⓓ</sup>, F. Bosco<sup>2 Ⓓ</sup>, C. Gieser<sup>3 Ⓓ</sup>, S. Suri<sup>4</sup>, J. C. Mottram<sup>2</sup>, R. Kuiper<sup>5 Ⓓ</sup>, T. Henning<sup>2</sup>, Á. Sánchez-Monge<sup>6,7,8,9 Ⓓ</sup>, H. Linz<sup>2 Ⓓ</sup>, R. E. Pudritz<sup>10 Ⓓ</sup>, D. Semenov<sup>2,11 Ⓓ</sup>, J. M. Winters<sup>12 Ⓓ</sup>, T. Möller<sup>9 Ⓓ</sup>, M. T. Beltrán<sup>13 Ⓓ</sup>, T. Csengeri<sup>14 Ⓓ</sup>, R. Galván-Madrid<sup>15 Ⓓ</sup>, K. G. Johnston<sup>16 Ⓓ</sup>, E. Keto<sup>17</sup>, P. D. Klaassen<sup>18 Ⓓ</sup>, S. Leurini<sup>19 Ⓓ</sup>, S. N. Longmore<sup>20 Ⓓ</sup>, S. L. Lumsden<sup>16</sup>, L. T. Maud<sup>21</sup>, L. Moscadelli<sup>13 Ⓓ</sup>, A. Palau<sup>15 Ⓓ</sup>, T. Peters<sup>22</sup>, S. E. Ragan<sup>23 Ⓓ</sup>, J. S. Urquhart<sup>24</sup>, Q. Zhang<sup>17 Ⓓ</sup>, and H. Zinnecker<sup>25</sup>

*(Affiliations can be found after the references)*

Received 29 November 2022 / Accepted 24 April 2023

## ABSTRACT

**Context.** The fragmentation mode of high-mass molecular clumps and the accretion processes that form the most massive stars ($M \ge 8 M_\odot$) are still not well understood. A growing number of case studies have found massive young stellar objects (MYSOs) to harbour disk-like structures, painting a picture that the formation of high-mass stars may proceed through disk accretion, similar to that of lower-mass stars. However, the properties of such structures have yet to be uniformly and systematically characterised.

**Aims.** The aim of this work is to uniformly study the kinematic properties of a large sample of MYSOs and characterise the stability of possible circumstellar disks against gravitational fragmentation.

**Methods.** We have undertaken a large observational programme (CORE) making use of interferometric observations from the Northern Extended Millimetre Array (NOEMA) for a sample of 20 luminous ($L > 10^4 L_\odot$) protostellar objects in the 1.37 mm wavelength regime in both continuum and spectral line emission, reaching 0.4" resolution (800 au at 2 kpc).

**Results.** We present the gas kinematics of the full sample and detect dense gas emission surrounding 15 regions within the CORE sample. Using the dense gas tracer CH$_3$CN, we find velocity gradients across 13 cores perpendicular to the directions of bipolar molecular outflows, making them excellent disk candidates. The extent of the CH$_3$CN emission tracing the disk candidates varies from 1800 to 8500 au. Analysing the free-fall to rotational timescales, we find that the sources are rotationally supported. The rotation profiles of some disk candidates are well described by differential rotation while for others the profiles are poorly resolved. Fitting the velocity profiles with a Keplerian model, we find protostellar masses in the range of ~10–25 $M_\odot$. Modelling the level population of CH$_3$CN (12$_K$–11$_K$) K = 0–6 lines, we present temperature maps and find median temperature in the range 70–210 K with a diversity in distributions. Radial profiles of the specific angular momentum (j) for the best disk candidates span a range of 1–2 orders of magnitude, on average ~$10^{-3}$ km s$^{-1}$ pc, and they follow $j \propto r^{1.7}$, which is consistent with a poorly resolved rotating and infalling envelope-disk model. Studying the Toomre stability of the disk candidates, we find almost all (11 out of 13) disk candidates to be prone to fragmentation due to gravitational instabilities at the scales probed by our observations, as a result of their high disk to stellar mass ratio. In particular, disks with masses greater than ~10–20% of the mass of their host (proto)stars are Toomre unstable, and more luminous YSOs tend to have disks that are more massive compared to their host star and hence more prone to fragmentation.

**Conclusions.** In this work, we show that most disk structures around high-mass YSOs are prone to disk fragmentation early in their formation due to their high disk to stellar mass ratio. This impacts the accretion evolution of high-mass protostars which will have significant implications for the formation of the most massive stars.

**Key words.** stars: formation – stars: massive – stars: kinematics and dynamics – techniques: interferometric

## 1. Introduction

There is a growing consensus that the formation of high-mass stars ($M \ge 8 M_\odot$) proceeds through disk accretion, similar to that of lower-mass stars. This has been shown through numerical simulations (e.g. Yorke & Sonnhalter 2002; Krumholz et al. 2009; Peters et al. 2010a; Kuiper et al. 2010, 2011; Kuiper & Yorke 2013; Klassen et al. 2016) and is supported by ubiquitous observations of molecular outflows in such regions (e.g. Henning

et al. 2000; Zhang et al. 2001, 2005; Beuther et al. 2002; Fallscheer et al. 2009; Leurini et al. 2011; Frank et al. 2014; Maud et al. 2015). Direct detection of disks around high-mass protostars has been difficult for many reasons. Firstly, high-mass star-forming regions are rare compared to their lower-mass counterparts (Kroupa 2001; Chabrier 2003), and they are typically located at further distances, requiring observations at a high angular resolution to study them in great detail. High-mass stars evolve much more rapidly than low-mass stars and, in their short formation phase, they are deeply embedded within their natal molecular cloud. Furthermore, a significant fraction of OB-type stars are found in binary or multiple systems (Apai et al. 2007; see the review by Offner et al. 2022). Spectroscopic studies have found >70% of main sequence O-type stars in close binary pairs,

<footnote>* Based on observations from an IRAM large programme L14AB. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).</footnote>

<footnote>** Data products for the CORE project can be found on the CORE collaboration website: https://www.mpia.de/core</footnote>
attributing this outcome to their formation process rather than
tidal capture or other scenarios (Chini et al. 2012; Sana et al.
2012). A pilot survey dedicated to the study of multiplicity in
massive young stellar objects (YSOs) using K-band observa-
tions has found a binary fraction of at least 50% (Pomohaci et al.
2019). High angular resolution observations are therefore espe-
cially critical for disentangling contributions from individual
protostars due to this clustered mode of formation of high-mass
stars (see review by Motte et al. 2018).

One of the first systematic searches for disks around
intermediate- and high-mass YSOs at high angular resolution
(tens of milli-arcsecond) was performed using the mid-infrared
instrument (MIDI) on the Very Large Telescope Interferom-
eter (VLTI) in the *N* band (8–13 µm; Boley et al. 2013).
Fitting the visibilities with 2D geometric models, Boley et al.
(2013) derived size, orientation, and elongation of the circum-
stellar material for a sub-sample of 20 objects. They found
that almost all objects showed deviations from spherical sym-
metry. In general, the compact mid-IR emission is elongated
along the direction of the circumstellar disk and perpendic-
ular to the outflow (e.g. S255 IRS3: Zinchenko et al. 2018).
However, sometimes the illumination footprints of outflow cones
can dominate the thermal emission at mid-IR, and therefore the
elongated nature of the intensity distribution at mid-IR for some
sources may be due to outflows rather than the embedded disks
(De Buizer et al. 2017), and hence the geometric models of the
MIDI data are susceptible to both disks and outflow cavity walls
in principle (Boley et al. 2013).

Since the majority of disk material is colder than what the
near- and mid-IR regimes can trace, (sub-)millimetre observa-
tions at high angular resolutions are necessary to obtain a full
picture, especially of the outermost regions of the disk and
to study the properties of the transition between the infalling
envelope and the disk. A review by Beltrán & de Wit (2016),
summarising the evidence for disks around intermediate- to
high-mass stars by compiling multi-wavelength observations
from the literature, confirmed that early B- to late O-type
angular resolution currently possible for 1.3 mm observations in the
northern sky (0.4", 800 au at 2 kpc) using the Northern Extended
Millimetre Array (NOEMA) interferometer and the IRAM 30-
m telescope (PI: Henrik Beuther). One of the objectives of this
large programme, CORE<fnref n="1" />, was to investigate the fragmentation
of cloud cores via gravitational Jeans instability (Beuther et al.
2018). We found a diversity in fragmentation of parsec-sized
clumps into cores with several thousand au sizes, whereby some
regions were found to have many cores, while other regions con-
tain a very low number of cores. We attributed this diversity
in fragmentation to a variety in initial density structures and/or
different initial magnetic field strengths.
The aim of the work presented in this paper is to study
the small-scale kinematics of the CORE sample to search for
and characterise disks around high-mass YSOs. In particular,
we study the stability of disk candidates against gravitational
instabilities in order to understand the importance of disk frag-
mentation in setting the final stellar mass distribution of a given
region. Other results from the CORE survey include various case
studies of individual regions (Ahmadi et al. 2018; Cesaroni et al.
2019; Bosco et al. 2019; Gieser et al. 2019; Mottram et al. 2020;
Moscadelli et al. 2021; Suri et al. 2021) as well as a chemical
analysis of the early phase of high-mass star formation based on
the full sample (Gieser et al. 2021).

It is important to target a large sample of high-mass YSOs
and uniformly search for and characterise the properties of a
statistically significant sample of disk candidates in order to
determine whether disk-mediated accretion is in fact the stan-
dard mechanism for high-mass star formation. In recent years,
a growing number of case studies have found disk-like struc-
tures around intermediate- to high-mass (proto)stars (Kraus et al.
2010; Johnston et al. 2015, 2020; Zapata et al. 2015, 2019; Ilee
et al. 2016, 2018; Chen et al. 2016; Cesaroni et al. 2017; Csengeri
et al. 2018; Ahmadi et al. 2018; Maud et al. 2018, 2019; Sanna
et al. 2019; Moscadelli et al. 2019; Guzmán et al. 2020; Tanaka
et al. 2020). If disk-mediated accretion is the standard mecha-
nism for high-mass star formation, then a closely related question
that needs investigation is whether or not the existence of close
binary or multiple systems is the product of disk fragmenta-
tion. Massive disks are prone to fragmentation via gravitational

instabilities due to high gas densities and accretion rates (Kratter & Matzner 2006; Kratter et al. 2010; review by Kratter & Lodato 2016). Recent numerical simulations have shown that binary and multiple systems with separations on the scales of hundreds of au can be formed as a result of disk fragmentation (Meyer et al. 2018; Oliva & Kuiper 2020), confirming speculations based on earlier, lower-resolution simulations that identified fragmentation as a crucial process in accretion flows around high-mass stars but that did not resolve the disk scale (Peters et al. 2010b). Observational studies have also started to resolve arc-, ring-, or spiral-like substructures that presumably develop as a result of gravitational instabilities in disks and inner envelopes (Maud et al. 2019; Zapata et al. 2019; Johnston et al. 2020). Furthermore, a recent study by Ramírez-Tannus et al. (2021) has shown that velocity dispersion of massive stars in young clusters increases as they age, suggesting that binaries are formed at larger separations and migrate inwards as a result of interactions with a remnant accretion disk or with other YSOs in the system. Therefore, both disk fragmentation and disk-mediated inward migration of fragments can be responsible for the observed close binary or multiple systems. Disk fragmentation may also be an important mechanism in seeding the initial binary and multiple systems before the inward migration occurs.

We have observed 20 high-mass YSOs at the highest angu-
lar resolution currently possible for 1.3 mm observations in the
northern sky (0.4", 800 au at 2 kpc) using the Northern Extended
Millimetre Array (NOEMA) interferometer and the IRAM 30-
m telescope (PI: Henrik Beuther). One of the objectives of this
large programme, CORE<fnref n="1" />, was to investigate the fragmentation
of cloud cores via gravitational Jeans instability (Beuther et al.
2018). We found a diversity in fragmentation of parsec-sized
clumps into cores with several thousand au sizes, whereby some
regions were found to have many cores, while other regions con-
tain a very low number of cores. We attributed this diversity
in fragmentation to a variety in initial density structures and/or
different initial magnetic field strengths.

The aim of the work presented in this paper is to study
the small-scale kinematics of the CORE sample to search for
and characterise disks around high-mass YSOs. In particular,
we study the stability of disk candidates against gravitational
instabilities in order to understand the importance of disk frag-
mentation in setting the final stellar mass distribution of a given
region. Other results from the CORE survey include various case
studies of individual regions (Ahmadi et al. 2018; Cesaroni et al.
2019; Bosco et al. 2019; Gieser et al. 2019; Mottram et al. 2020;
Moscadelli et al. 2021; Suri et al. 2021) as well as a chemical
analysis of the early phase of high-mass star formation based on
the full sample (Gieser et al. 2021).

The paper is organised as follows. In Sect. 2 we describe how
the CORE sample was selected, followed by an overview of the
observations and data reduction techniques in Sect. 3. We out-
line the observational results in Sect. 4 with a brief description
of the large-scale core fragmentation and kinematics, and then
we turn our focus to the dense gas distribution and kinematics.
In Sect. 5 we begin with an analysis of the temperature distribu-
tion and provide estimates for gas masses. We continue with a
study of the candidacy of the sources as disks by studying their
rotation profiles, timescales, and specific angular momenta, and
finally investigate the stability of the disk candidates. The paper
concludes with a summary of our findings in Sect. 6.

<footnote n="1">¹ http://www.mpia.de/core</footnote>
<table><caption>**Table 1.** Positions and properties of the CORE sample, grouped in track-sharing pairs, adapted from Beuther et al. (2018).</caption><thead><tr><th>Source</th><th>R.A.<br>(J2000.0)</th><th>Dec.<br>(J2000.0)</th><th>Distance<br>(kpc)</th><th>Luminosity<br>(10<sup>4</sup> L<sub>☉</sub>)</th><th>σ<sub>cont</sub><sup>(a)</sup><br>(mJy beam<sup>-1</sup>)</th><th>IR-bright</th><th>a.f.<sup>(b)</sup></th><th>Ref.<sup>(c)</sup></th></tr></thead><tbody><tr><td>IRAS 23151+5912</td><td>23:17:21.01</td><td>+59:28:47.49</td><td>3.3</td><td>2.4</td><td>0.19</td><td>+</td><td>w</td><td>d1<sup>(*)</sup>, l2</td></tr><tr><td>IRAS 23033+5951</td><td>23:05:25.00</td><td>+60:08:15.49</td><td>4.3</td><td>1.7</td><td>0.46</td><td>—</td><td>cm, w</td><td>d2, l1</td></tr><tr><td>AFGL2591</td><td>20:29:24.86</td><td>+40:11:19.40</td><td>3.3</td><td>20.0</td><td>0.60</td><td>+</td><td>cm, w</td><td>d3<sup>(*,†)</sup>, l1</td></tr><tr><td>G75.78+0.34</td><td>20:21:44.03</td><td>+37:26:37.70</td><td>3.8</td><td>11.0</td><td>0.60</td><td>—</td><td>cm, m</td><td>d4<sup>(*,†)</sup>, l1</td></tr><tr><td>S87 IRS1</td><td>19:46:20.14</td><td>+24:35:29.00</td><td>2.2</td><td>2.5</td><td>0.23</td><td>+</td><td>cm</td><td>d5<sup>(†)</sup>, l1</td></tr><tr><td>S106</td><td>20:27:26.77</td><td>+37:22:47.70</td><td>1.3</td><td>3.4</td><td>1.25</td><td>+</td><td>cm, w</td><td>d6<sup>(*)</sup>, l2</td></tr><tr><td>IRAS 21078+5211</td><td>21:09:21.64</td><td>+52:22:37.50</td><td>1.5</td><td>1.3</td><td>0.60</td><td>—</td><td>cm, w</td><td>d11<sup>(*)</sup></td></tr><tr><td>G100.3779-03.578</td><td>22:16:10.35</td><td>+52:21:34.70</td><td>3.5</td><td>1.5</td><td>0.08</td><td>+</td><td>w</td><td>d1<sup>(*)</sup>, l2</td></tr><tr><td>G084.9505-00.691</td><td>20:55:32.47</td><td>+44:06:10.10</td><td>5.5</td><td>1.3</td><td>0.10</td><td>+</td><td>w</td><td>d2, l2</td></tr><tr><td>G094.6028-01.797</td><td>21:39:58.25</td><td>+50:14:20.90</td><td>4.0</td><td>2.8</td><td>0.14</td><td>+</td><td>w, m</td><td>d7<sup>(*)</sup>, l2</td></tr><tr><td>Cep A HW2</td><td>22:56:17.98</td><td>+62:01:49.50</td><td>0.8</td><td>1.5</td><td>4.00</td><td>—</td><td>cm, w, m</td><td>d8<sup>(†)</sup>, l1</td></tr><tr><td>NGC 7538 IRS9</td><td>23:14:01.68</td><td>+61:27:19.10</td><td>2.7</td><td>2.3</td><td>0.30</td><td>+</td><td>w</td><td>d9<sup>(*,†)</sup>, l1</td></tr><tr><td>W3(H<sub>2</sub>O)</td><td>02:27:04.60</td><td>+61:52:24.73</td><td>2.0</td><td>8.3</td><td>4.50</td><td>—</td><td>cm, w, m</td><td>d10<sup>(*)</sup>, l2</td></tr><tr><td>W3 IRS4</td><td>02:25:31.22</td><td>+62:06:21.00</td><td>2.0</td><td>4.5</td><td>0.60</td><td>+</td><td>cm, w</td><td>d10, l1</td></tr><tr><td>G108.7575-00.986</td><td>22:58:47.25</td><td>+58:45:01.60</td><td>4.3</td><td>1.4</td><td>0.25</td><td>+</td><td>w, m</td><td>d2, l3</td></tr><tr><td>IRAS 23385+6053</td><td>23:40:54.40</td><td>+61:10:28.02</td><td>4.9</td><td>1.6</td><td>0.25</td><td>—</td><td>w</td><td>d12</td></tr><tr><td>G138.2957+01.555</td><td>03:01:31.32</td><td>+60:29:13.20</td><td>2.9</td><td>1.4</td><td>0.16</td><td>+</td><td>cm, w</td><td>d2<sup>(†)</sup>, l1</td></tr><tr><td>G139.9091+00.197</td><td>03:07:24.52</td><td>+58:30:48.30</td><td>3.2<sup>(‡)</sup></td><td>1.1</td><td>0.17</td><td>+</td><td>cm, w</td><td>d2, l1</td></tr><tr><td colspan="9"><strong>Pilot study</strong></td></tr><tr><td>NGC 7538 IRS1</td><td>23:13:45.36</td><td>+61:28:10.55</td><td>2.7</td><td>21.0</td><td>10.0</td><td>+</td><td>cm, w, m</td><td>d9<sup>(*,†)</sup>, l1</td></tr><tr><td>NGC 7538 S</td><td>23:13:44.86</td><td>+61:26:48.10</td><td>2.7</td><td>1.5</td><td>0.60</td><td>—</td><td>w, m</td><td>d9<sup>(*,†)</sup>, l1</td></tr></tbody></table>

**Notes.** (a) rms noise in the continuum maps imaged with a robust 0.1 (uniform) weighting from Beuther et al. (2018). (b) Associated features (a.f.): cm: cm continuum; w: H<sub>2</sub>O maser; m: CH<sub>3</sub>OH maser. (c) References for distances and luminosities: d1: Choi et al. (2014), d2: Urquhart et al. (2011), d3: Rygl et al. (2012), d4: Ando et al. (2011), d5: Xu et al. (2009), d6: Xu et al. (2013), d7: Choi et al. (2014), Sakai et al. (2019), d8: Xu et al. (2016), d9: Moscadelli et al. (2009), d10: Hachisuka et al. (2006), Xu et al. (2006), d11: Molinari et al. (1996), d12: Molinari et al. (1998); l1: rms survey database (http://rms.leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi) using SED fitting from Mottram et al. (2011) including Herschel fluxes and the latest distance determination; l2: rms survey database, using SED fitting from Mottram et al. (2011) updated to the latest distance determination; l3: rms survey database, calculated from the MSX 21 µm flux using the scaling relation derived by Mottram et al. (2011) and updated to the latest distance determination. (*)Distances in agreement with Reid et al. (2019). (†)Distances in agreement with Mège et al. (2021). (‡)Mège et al. (2021) report a kinematic distance of 1.6 ± 0.7 kpc.

## 2. Source selection

As outlined in Beuther et al. (2018) in more detail, the CORE sample of 20 sources includes all targets within the Red MSX (RMS) survey<fnref n="2" /> (Lumsden et al. 2013) for which high-resolution millimetre data did not previously exist and that conform to the following criteria. To ensure that the targets are high-mass star-forming regions, we made a luminosity cut at $L > 10^4 L_☉$, which corresponds to a star with a mass of at least 8 $M_☉$. To be able to observe the sample at high resolution with good image fidelity, a distance cut at $d < 6$ kpc was enforced and only sources with declination $>24^\circ$ were included to ensure good uv-coverage. To target the early phase of formation with active accretion, we only selected sources whose fluxes at 21 µm are at least two times greater than their 8 µm flux, using the RMS survey, since it is expected that the IR colour evolves from red<fnref n="3" /> to blue. In this way, the final sample size was narrowed down to 18 regions along with two additional sources (NGC 7538 IRS1 and NGC 7538 S)

which were observed previously with almost the same setup; we subsequently refer to these two sources as the CORE pilot regions (Beuther et al. 2012; 2013; Feng et al. 2016). Table 1 lists the names and positions of the sources with their properties, including features associated with the regions such as the presence of various masers and/or centimetre continuum emission as proxy for the star formation activity. The threshold between IR-bright and IR-dark is defined according to Motte et al. (2007) based on cores in the Cygnus X star-forming region at a distance $d = 1.7$ kpc, such that if a core with a luminosity of 1000 $L_☉$ has 21 µm flux ($S_{21µm}$) below 10 Jy it would be categorised as IR quiet. Therefore, regions with $S_{21µm} < 10$ Jy ($\frac{1.7 \text{ kpc}}{d}^2$) ($\frac{L}{1000 L_☉}$) are classified as IR-dark.

## 3. Observations and data reduction

Observations in the 1.3 mm wavelength regime of the CORE project began in June 2014 and finished in January 2017, consisting of a total of more than 400 hours of observations with originally the Plateau de Bure Interferometer (PdBI) that evolved

<footnote n="2">² http://rms.leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi</footnote>

<footnote n="3">³ When the flux at 21 µm is greater than at 8 µm.</footnote>
<table><caption>**Table 2.** Frequency setup of the narrow-band correlator and important lines covered.</caption><thead><tr><th>Spectral unit</th><th>Frequency range (MHz)</th><th>Important lines</th></tr></thead><tbody><tr><td>L01</td><td>220 690.6–220 769.7</td><td>CH<sub>3</sub>CN (12<sub>K</sub>–11<sub>K</sub>) K = 0–3</td></tr><tr><td>L02</td><td>220 630.6–220 709.7</td><td>CH<sub>3</sub>CN (12<sub>K</sub>–11<sub>K</sub>) K = 4–5, C<sub>2</sub>H<sub>5</sub>CN (25<sub>22,4</sub> – 24<sub>22,3</sub>)</td></tr><tr><td>L03</td><td>220 570.6–220 649.7</td><td>CH<sub>3</sub>CN (12<sub>K</sub>–11<sub>K</sub>) K = 5–6, CH<sub>3</sub><sup>13</sup>CN (12<sub>K</sub>–11<sub>K</sub>) K = 0–3, HNCO (10<sub>1,9</sub>–9<sub>1,8</sub>)</td></tr><tr><td>L04</td><td>220 130.6–220 209.7</td><td>HCOOCH<sub>3</sub> (17<sub>4,13</sub>–16<sub>4,12</sub>) A/E, H<sub>2</sub>CCO (11<sub>1,11</sub>–10<sub>1,10</sub>)</td></tr><tr><td>L05</td><td>218 860.6–218 939.7</td><td>OCS (18 – 17)</td></tr><tr><td>L06</td><td>218 415.6–218 494.7</td><td>H<sub>2</sub>CO (3<sub>2,2</sub>–2<sub>2,1</sub>), CH<sub>3</sub>OH (4<sub>2</sub>–3<sub>1</sub>)</td></tr><tr><td>L07</td><td>218 280.6–218 359.7</td><td>HC<sub>3</sub>N (24–23), HCOOCH<sub>3</sub> (17<sub>3,14</sub>–16<sub>3,13</sub>) A</td></tr><tr><td>L08</td><td>218 180.6–218 259.7</td><td>H<sub>2</sub>CO (3<sub>0,3</sub>–2<sub>0,2</sub>), O<sup>13</sup>CS (18–17)</td></tr></tbody></table>

into NOEMA over the time of the project. We observed in three different configurations of the interferometer: A-configuration being the most extended, D-configuration being the most compact array, and B-configuration covering baselines between the two. Baselines roughly in the range of 15–765 m were covered, therefore the NOEMA observations are not sensitive to structures larger than ~15" (0.1 pc at 2 kpc) at 220 GHz. The exact baseline ranges for all tracks are listed in Table 3 of Beuther et al. (2018). On-source observations were taken in roughly 20-min increments distributed over an observing run (i.e. a track) and interleaved with observations of various calibration sources. We observed in track-sharing mode meaning two science targets were observed in one track, with one track corresponding to a total of roughly 8 h of observations. We paired sources with similar positions in the sky together in order to be able to use the same calibrators. With the sources being observed in three different configurations, in total at least three half tracks were observed per source. In many cases more tracks were needed to achieve the required sensitivity, better uv-coverage, and depending on the weather conditions. The observations of the two pilot sources were carried out in early 2011 in the A- and B-array configurations of NOEMA only (Beuther et al. 2012).

The spectral setup of the CORE survey was outlined in Ahmadi et al. (2018). To summarise, we observed the frequency range of 217 078.6–220 859.5 MHz with the wide-band correlator WideX, in both horizontal and vertical polarisations with a fixed spectral resolution of 1.95 MHz (~2.7 km s⁻¹ at 219 GHz). For higher spectral resolution, we used 8 narrow-band correlator units with 80 MHz bandwidth and spectral resolution of 0.312 MHz (~0.43 km s⁻¹). The exact bandwidth coverage of this correlator and the important lines observed are summarised in Table 2. The main species used in this work were observed with the narrow-band correlator, smoothed to 0.5 km s⁻¹ spectral resolution. Even though the observations were carried out over a period of a few years when the NOEMA interferometer was upgraded from 6 to 8 antennas, the narrow-band correlator was only able to accept the signal from the 6 antennas that provided the best uv-coverage.

Calibration of the NOEMA observations was performed with the CLIC programme of the Grenoble Image and Line Data Analysis Software (GILDAS<fnref n="4" />) package developed by IRAM and Observatoire de Grenoble. In most cases, two phase and amplitude calibrators were observed. Strong quasars (e.g. 3C84, 3C273, or 3C454.3) were used for the bandpass calibration. For most of the target observations, the source MWC349 was used for absolute flux calibration, assuming a model flux of 1.86 Jy

at 220 GHz. For a few tracks, the source LKHα 101 was used for the flux calibration. The absolute flux scale is correct to within 20%, according to the observatory standard.

The imaging of the NOEMA observations was performed with the MAPPING programme of GILDAS. The continuum was extracted from line-free regions of the WideX band as described in Beuther et al. (2018) and was imaged with uniform weighting yielding synthesised beam sizes of ~0.32"–0.5". The rms noise of the continuum images are listed in Table 1 and depend on the prominence of sidelobes in the images whereby the brighter sources tend to be noisier. To study the detailed kinematics of the regions on small scales, we analysed the highest spectral resolution molecular line observations from the narrow-band correlator listed in Table 2. Each spectral cube was resampled to have 60 channels centred at the frequency of the given transition with a spectral resolution of 0.5 km s⁻¹. For each molecular line transition, the continuum was subtracted in the uv-plane by selecting the line-free channels in the spectral unit in which the transition lies. The only exceptions were molecular transitions in spectral unit L03 because too few line-free channels were present. For this spectral unit, we used line-free channels from spectral unit L02 assuming there is no significant spectral slope across these adjacent units. The resulting visibilities were then imaged with the CLARK algorithm (Clark 1980) with different weightings, discussed further in the analysis (see Sect. 4.2).

Observations with the IRAM 30-m telescope were taken in March and May of 2015 using the EMIR multi-band mm-wave receiver (Carter et al. 2012). On-the-fly maps (whereby the telescope continuously slews smoothly and rapidly across a field) covering an area of 1.5' by 1.5' were observed for all 20 sources centred at the phase centre of the NOEMA observations. The half-power-beam-width of the observations is 11.3" at 217 GHz and 10.6" at 231 GHz. Two 4 GHz wide bands were positioned at 215, and 219 GHz in the lower sideband with a small overlap to cover the spectral setup of the NOEMA observations with the WIDEX correlator. Another two 4 GHz wide units on the upper sideband were centred at 231 and 234 GHz, also with a small overlap. A Fourier Transform Spectrometer (FTS) backend was employed with a channel width of 195 kHz corresponding to a resolution of 0.27 km s⁻¹ at 217 GHz and 0.25 km s⁻¹ at 233 GHz. The single-dish data were reduced using the CLASS programme of GILDAS and the procedure is detailed in Appendix A.1 of Mottram et al. (2020). In the work presented here, the single-dish observations were only used to infer outflow directions to strengthen the existence of disks using typical outflow and shock tracers. For this purpose, we also merged the single-dish and NOEMA observations for

<footnote n="4">⁴ http://www.iram.fr/IRAMFR/GILDAS</footnote>
$^{13}$CO(2–1) using the MAPPING software and produced images with natural weighting, and the Steer-Dewdney-Ito (SDI; Steer et al. 1984) algorithm which is more appropriate when the distribution of emission is extended. The spectral resolution of this merged cube is 3.0 km s$^{-1}$. More details on the merging process can be found in Appendix B of Mottram et al. (2020) with an extensive exploration of the effects of using different imaging algorithms and weightings.

# 4. Observational results

## 4.1. Core fragmentation, kinematics, and evolutionary stage

We detect strong millimetre continuum emission from all sources within the CORE survey, mainly associated with cold dust. Contours in Fig. 1 show the continuum observations where a diversity of fragmentation levels among the sources can be seen, with some regions harbouring as many as 20 cores (e.g. IRAS 21078, while others remain singular (e.g. AFGL2591). To understand this diversity, we performed a minimum-spanning-tree analysis and found the separation between the cores to be of the order of the thermal Jeans length or smaller (Beuther et al. 2018). This analysis indicated that the initial density structure or magnetic fields may be among the important factors that set the fragmentation diversity we see on these clump scales. Because the smallest separation between the cores were often found to be at our resolution limit, it is expected that further fragmentation can be found on even smaller scales, likely as a result of disk fragmentation which is the focus of this paper.

We made use of H$_2$CO and CH$_3$OH transitions that are excited at lower gas temperatures and densities than CH$_3$CN to probe the distribution of gas on large scales. Figure 1 shows the intensity-weighted mean velocity (first moment) maps of H$_2$CO (3$_{0.3}$–2$_{0.2}$) for the entire sample. We detect gas with velocities varying of the order of a few km s$^{-1}$ in each field. For W3IRS4 that has a filamentary structure, we find a velocity gradient along the filament, consistent with flows of material along the filament onto the hot core that resides at the end of the filament (Mottram et al. 2020). There exist a separate velocity gradient across the hot core with a different position angle and amplitude than the larger scale flow along the filament. This finding is consistent with the filamentary paradigm of star formation (see reviews by André et al. 2014 and Hacar et al. 2022).

We further find a variety of stages in the evolutionary sequence of high-mass star formation within some of the targeted fields, hinting at the importance of the environment on triggering star formation. For example, in the case of W3IRS4 mentioned above, we find two cold cores along a filamentary structure and a hot core situated at the junction of the filament and an expanding UCH II region known as W3C. The UCH II region is seen as a circular structure to the south-west with almost no H$_2$CO emission in the interferometric observations. We can compare this source to another UCH II region in the CORE sample, W3(OH), which is at the same distance and has a similar mass sensitivity to W3IRS4. W3(OH) is the western-most feature seen in the map labelled as W3(H$_2$O), around which we detect extended H$_2$CO emission, unlike for W3C. This difference can in part be attributed to the fact that W3C is probably at a later evolutionary stage than W3(OH) as suggested by its larger extent, hence photoionisation may have already destroyed many of the species that would often be found during an earlier evolutionary phase.

Although we attempted to observe regions at the same evolutionary stage, the CORE sample does inevitably span a range in evolutionary stage, as high-mass star formation occurs on much

shorter timescales compared to the formation of low-mass stars. This is evident because some YSOs show a rich hot core chemistry (e.g. AFGL2591, G75.78, Cep A HW2, and W3(H$_2$O)), while others are weaker in line emission (e.g. S87 IRS1, S106, G108, G138, G139). For two of these weak sources, namely S106 and G139, we barely detect any H$_2$CO (3$_{0.3}$–2$_{0.2}$) emission at the position of the cores (see Fig. 1). Furthermore, there is no correlation between the richness in lines and the infrared (IR) brightness of the sources (see Table 1). The fact that some of the most chemically rich sources in our survey are IR-dark points to their deeply embedded nature, during which they are already capable of producing a rich chemistry and driving powerful outflows. We refer to Gieser et al. (2021) for a more comprehensive study of the observed chemistry in the CORE sample.

## 4.2. Dense gas kinematics

To study candidate disks, one needs the right spectral line tracer that is chemically abundant, excited only in the innermost regions of the core, and not optically thick. CH$_3$CN has proven to be an excellent dense gas and disk tracer, and using its higher excited transitions along with its rarer CH$_3$$^{13}$CN isotopologue, line optical thickness issues can be circumvented (e.g. Sánchez-Monge et al. 2014; Cesaroni et al. 2017; Beuther et al. 2018). Additionally, due to its molecular structure, CH$_3$CN splits into a ladder of closely spaced transitions (i.e. K-ladder), the level populations of which can provide a measure for core temperatures and densities. In this work, we make use of CH$_3$CN (12$_K$–11$_K$) K = 0–6 transitions with upper energies in the range of 69 to 325 K to study the dense gas kinematics and physical conditions of the full sample. Initially, we imaged the CH$_3$CN cubes for the entire sample with a uniform weighting (robust 0.1 in GILDAS) to reach the highest possible angular resolution. For sources for which CH$_3$CN was either not detected or only one or two of the lowest K transitions were detected, we reimaged the data cubes with weightings that help in recovering more extended emission but that result in coarser spatial resolution (a higher robust parameter closer to natural weighting). The sources for which up-weighting the compact configuration helped in recovering some CH$_3$CN emission were IRAS 23033, IRAS 23385, G100, and G084. Even with such a strategy, we could not detect CH$_3$CN emission for 5 of the 20 core sources: S87 IRS1, S106, G108, G138, and G139. Interestingly, all of these objects are IR-bright confirming our hypothesis that these sources are likely at a slightly more advanced evolutionary phase in which the warmer environment has already destroyed this species. It is still possible that these sources may harbour disks that are at a later evolutionary stage that cannot be traced by CH$_3$CN. These five sources were therefore omitted from our search for circumstellar disks.

In Fig. 2 we show the intensity-weighted mean velocity (first moment) maps of CH$_3$CN (12$_K$–11$_K$) for the remaining 15 sources in the CORE sample in search of velocity gradients that trace rotation. The details of image weighting parameters and the associated beam sizes, as well as the ranges over which the emission is integrated are listed in Table 3. The corresponding zeroth and second moment maps can be found in Appendix A. The threshold used for creating the moment maps is set to $6\sigma$, where $\sigma$ is estimated as the root mean square noise in an emission-free area in the channel corresponding to the peak of the line; this is also listed in Table 3. We detect velocity gradients across almost all cores with magnitudes of a few km s$^{-1}$, depicted by black dotted lines in the panels, spanning a physical projected distance of 2000–8000 au. The bulk of the
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Intensity-weighted mean velocity (first moment) maps of H₂CO (3₀.₃ – 2₀.₂) showing the large-scale kinematics of the full CORE sample. The contours correspond to the continuum maps imaged with uniform weighting as presented in Beuther et al. (2018). The outermost three contours correspond to 5, 10, and 20σ levels, then increasing in steps of 15σ. For IRAS 23151, NGC 7538 IRS1, NGC 7538 IRS9 and AFGL2591 the outermost three contours correspond to 5, 15, and 40σ levels, then increasing in steps of 25σ (see Table 3 for σ values). The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel. When displayed in Adobe Acrobat, it is possible to switch to the [first moment map of CH₃OH(4'–3')] and back to the [first moment map of H₂CO (3₀.₃ – 2₀.₂)] for comparison.</figcaption></figure>
<figure><img src="image_3.png" /></figure>

**Fig. 2.** Intensity-weighted mean velocity (first moment) maps of CH₃CN (12₃⁻11₃) showing the dense gas kinematics for 15 of the 20 sources in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively (see Sect. 4.3). The dotted lines indicate the position of the strongest velocity gradient tracing the disk, i.e. perpendicular to the rotation axis. The edges of the assumed disk extent are marked with an ×. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel. When displayed in Adobe Acrobat, it is possible to switch to the [decomposed velocity map] and back to the [first moment map] for comparison; otherwise, see Fig. B.2.
<table><caption>**Table 3.** Observational parameters for the disk candidates within the CORE survey.</caption><thead><tr><th>Source</th><th>Weight<sup>(a)</sup></th><th>Beam<br>("x", PA)</th><th>Res.<sup>(b)</sup><br>(au)</th><th>σ<sub>line</sub><sup>(c)</sup><br>(mJy beam<sup>-1</sup>)</th><th>Outflow<br>PA<sup>(d)</sup><br>(°)</th><th>Disk<br>PA<sup>(d)</sup><br>(°)</th><th>R<sub>blue</sub><sup>(e)</sup><br>(au)</th><th>R<sub>red</sub><sup>(f)</sup><br>(au)</th><th>Extent<sup>(g)</sup><br>(au)</th><th>v<sub>LSR</sub><sup>(h)</sup><br>(km s<sup>-1</sup>)</th><th>Integration range<sup>(i)</sup><br>(km s<sup>-1</sup>)</th><th><sup>12</sup>C/<sup>13</sup>C<sup>(j)</sup></th></tr></thead><tbody><tr><td>IRAS 23151</td><td>robust 0.1</td><td>0.42 × 0.35, 47°</td><td>1280</td><td>3.60</td><td>230</td><td>140</td><td>2680</td><td>2680</td><td>5360</td><td>-53.1</td><td>-64 → -48</td><td>93</td></tr><tr><td>IRAS 23033 MMS1b</td><td>robust 1</td><td>0.49 × 0.43, 9°</td><td>1980</td><td>3.86</td><td>140</td><td>50</td><td>520</td><td>1300</td><td>1820</td><td>-54.4</td><td>-59 → -48</td><td>81</td></tr><tr><td>AFGL2591</td><td>robust 0.1</td><td>0.45 × 0.34, -115°</td><td>1310</td><td>5.91</td><td>100</td><td>45</td><td>2680</td><td>3400</td><td>6080</td><td>-5.6</td><td>-11 → 0</td><td>68</td></tr><tr><td>G75.78</td><td>robust 0.1</td><td>0.45 × 0.34, 59°</td><td>1510</td><td>5.19</td><td>35</td><td>130</td><td>3920</td><td>3920</td><td>7840</td><td>0.2</td><td>-7 → 5</td><td>67</td></tr><tr><td>IRAS 21078</td><td>robust 0.1</td><td>0.45 × 0.31, 39°</td><td>570</td><td>4.04</td><td>45</td><td>100</td><td>1220</td><td>990</td><td>2210</td><td>-7.0</td><td>-15 → 2</td><td>69</td></tr><tr><td>G094</td><td>robust 0.1</td><td>0.38 × 0.36, -48°</td><td>1480</td><td>3.92</td><td>60</td><td>170</td><td>2900</td><td>2540</td><td>5440</td><td>-42.9</td><td>-52 → -35</td><td>80</td></tr><tr><td>Cep A</td><td>robust 0.1</td><td>0.40 × 0.36, 77°</td><td>300</td><td>6.94</td><td>50</td><td>110</td><td>930</td><td>930</td><td>1860</td><td>-5.6</td><td>-17 → 2</td><td>70</td></tr><tr><td>NGC 7538 IRS9</td><td>robust 0.1</td><td>0.40 × 0.35, 78°</td><td>1020</td><td>4.27</td><td>90</td><td>0</td><td>1725</td><td>770</td><td>2495</td><td>-56.9</td><td>-68 → -47</td><td>77</td></tr><tr><td>W3(H<sub>2</sub>O) E</td><td>robust 0.1</td><td>0.42 × 0.30, 87°</td><td>720</td><td>6.81</td><td>60</td><td>130</td><td>1930</td><td>2970</td><td>4900</td><td>-51.6</td><td>-57 → -39</td><td>79</td></tr><tr><td>W3(H<sub>2</sub>O) W</td><td>robust 0.1</td><td>0.42 × 0.30, 87°</td><td>720</td><td>6.81</td><td>30</td><td>100</td><td>1930</td><td>2520</td><td>4450</td><td>-48.7</td><td>-57 → -39</td><td>79</td></tr><tr><td>IRAS 23385</td><td>robust 5</td><td>0.93 × 0.88, 113°</td><td>4440</td><td>5.02</td><td>10</td><td>100</td><td>3980</td><td>4700</td><td>8680</td><td>-50.6</td><td>-57 → -43</td><td>90</td></tr><tr><td colspan="14"><strong>Pilot study</strong></td></tr><tr><td>NGC 7538 IRS1</td><td>robust 0.1</td><td>0.41 × 0.36, 118°</td><td>1030</td><td>3.88</td><td>130</td><td>30</td><td>2980</td><td>1990</td><td>4970</td><td>-56.3</td><td>-66 → -45</td><td>78</td></tr><tr><td>NGC 7538 S</td><td>robust 0.1</td><td>0.31 × 0.29, 104°</td><td>800</td><td>3.12</td><td>135</td><td>175</td><td>990</td><td>2190</td><td>3180</td><td>-61.4</td><td>-66 → -50</td><td>78</td></tr></tbody></table>

**Notes.** Images are from a combination of configurations A, B, and D arrays of NOEMA with the exception of pilot sources that lack D-array observations. <sup>(a)</sup>Robust weighting parameter using GILDAS definition whereby robust 0.1 corresponds to uniform weighting and robust 10 is closest to natural weighting. <sup>(b)</sup>Average linear resolution, using distances from Table 1. <sup>(c)</sup>rms noise in the emission-free region in the channel that has the peak of emission for CH<sub>3</sub>CN (12<sub>4</sub>–11<sub>4</sub>), therefore the unit is per 0.5 km s<sup>-1</sup> channel. <sup>(d)</sup>Measured counterclockwise from the north-south axis. <sup>(e)</sup>Extent of the blueshifted emission for CH<sub>3</sub>CN (12<sub>3</sub>–11<sub>3</sub>). <sup>(f)</sup>Extent of the redshifted emission for CH<sub>3</sub>CN (12<sub>3</sub>–11<sub>3</sub>). <sup>(g)</sup>Sum of R<sub>blue</sub> and R<sub>red</sub>. <sup>(h)</sup>Peak velocity of CH<sub>3</sub>CN (12<sub>3</sub>–11<sub>3</sub>) at the position of the continuum peak. <sup>(i)</sup>Integration ranges for the moment maps shown in Fig. 2. <sup>(j)</sup>Isotopic ratios calculated according to Wilson & Rood (1994), (see text for details).
emission for most sources is covered by many beams, allowing us to properly resolve these objects, with the exception of the sources for which we degraded the resolution in order to detect CH3CN. In these cases, the extent of the emission is only one or two beams across. For two of these sources in particular (G084 and G100), the velocity gradients in the first moment maps are not particularly convincing. Furthermore, IRAS 23033 consists of two cores (MMS1a and MMS1b), both of which have detection of CH3CN, with MMS1b the brighter of the two. Similarly, W3(H2O) consists of two cores, E and W corresponding to the eastern and western cores, respectively, approximately 2000 au apart (Ahmadi et al. 2018). Therefore, we detect CH3CN for 17 cores in total in this early stage of star formation.

Intensity-weighted first moment maps are inherently sensitive to noise in the spectrum and any asymmetries in the line profile (see Teague & Foreman-Mackey 2018 for a short overview). In order to overcome this issue and detect the velocity field beyond the 6σ threshold used in the creation of moment maps, which can often mask more emission than needed, we have fitted Gaussian profiles to the spectra of all sources using the Spectral-Cube tool within the PySpecKit package in Python (Ginsburg & Mirocha 2011; Ginsburg et al. 2019, 2022). Fitting one Gaussian profile assumes that there is only one velocity component at each position. In these so-called decomposed velocity maps, we detect emission that is slightly more extended in all directions for all sources (see Fig. B.2). For G084 and G100, velocity gradients become slightly clearer than in the first moment maps, albeit not as strong as the other sources, with amplitudes below 1 km s⁻¹, making them less reliable disk candidates at the resolution of our observations. Furthermore, the velocity gradient across IRAS 23033 MMS1a is in the same direction as the molecular outflow, depicted by blue and red arrows in Fig. 2 and discussed further in the next section (see CORE case study by Bosco et al. (2019) for a detailed analysis of IRAS 23033). Moreover, Mottram et al. (2020) performed an in-depth study of W3IRS4 using the CORE data and detected velocity gradients in a number of species around this core, but found that the axis of the velocity gradient changes depending on the species, hence likely tracing a different combination of infall, rotation, and outflow, depending on the tracer. They further concluded that any disk around W3IRS4 must be smaller than the resolution of our observations for this source (600 au). For these reasons, we remove G084, G100, IRAS 23033 MMS1a, and W3IRS4 from our list of disk candidates. For the remaining 13 best disk candidates, the velocity gradients observed using the dense gas tracer CH3CN are consistent and coherent with the motion of more diffuse gas on larger scales, for example as traced by H2CO (see Fig. 1). The position angles of the CH3CN velocity gradients are listed in Table 3, defined with respect to the north-south axis and measured counterclockwise.

### 4.3. Molecular outflows

We make use of outflow tracing molecules, mainly the *J* = 2–1 transition of CO (from our single-dish IRAM 30-m observations) and its ¹³CO isotopologue (from our single-dish, NOEMA, and merged observations), to infer the directions of molecular outflows being driven by these YSOs. In Appendix C we show the intensity maps of the blueshifted and redshifted gas with respect to the systemic velocity of the source by integrating over the blue and red wings of the emission for these transitions as well as for C¹⁸O (2–1) and SO (6$_5$–5$_4$). Since the interferometer filters out most of such extended emission, we make use of

the IRAM 30-m data with a field of view of 1.5' by 1.5' (and half-power beam width of ~11'' at this frequency) to infer the structure of these large-scale flows. Furthermore, the interferometer and in particular the combination of the interferometric with single-dish observations for ¹³CO (2–1) are useful in disentangling the ejection points of the outflows, especially in cores which have substructure as multiple outflows may exist, being driven by other objects (e.g. the case of W3(H₂O) presented in Ahmadi et al. 2018).

We observe extended high-velocity CO (2–1) and ¹³CO (2–1) emission in the vicinity of the most luminous source in the field for all 17 cores where we detect CH3CN emission (see Appendix C). In many of the sources, we detect clear bipolar molecular outflows (e.g. IRAS 23151, AFGL2591, Cep A, W3(H₂O)) while the case is more complex for others. As an example, for NGC 7538 IRS1, our data suggest a bipolar outflow in the NW-SE direction confirmed by other works (e.g. Keto 1991; Davis et al. 1998; Klaassen et al. 2011; Qiu et al. 2011; Beuther et al. 2013), while an ionised jet in the N-S direction has been detected by others (e.g. Gaume et al. 1995; Sandell et al. 2009; Moscadelli & Goddi 2014; Goddi et al. 2015). This ambiguity can be attributed to the fact that further fragmentation has been reported within this source by higher resolution observations in the centimetre regime (Goddi et al. 2015; Beuther et al. 2017), and made more complex because the central source is already ionising a hypercompact H II region. In fact, NGC 7538 IRS1 is the only source within the CORE survey where we detect the molecular line transitions in absorption against the strong continuum background at the peak position of the continuum emission spanning an area roughly the size of a synthesised beam (Beuther et al. 2012). Furthermore, for sources that are in complex environments such as IRAS 23033 and W3IRS4, establishing an exact direction for the molecular outflow(s) is difficult. For IRAS 23033, we adopt the directions of outflows determined by Bosco et al. (2019) who performed a detailed analysis of this region using the same CORE data. Similarly, for W3IRS4, Mottram et al. (2020) provided an in-depth study of this region using the CORE data and adopted an outflow direction roughly in the NNE-SSW direction, which is what we adopt here as well.

The estimated directions of blueshifted and redshifted emission of molecular outflows are shown on Fig. 2 as blue and red arrows, respectively. The outflow position angles are also summarised in Table 3. It is important to note that due to the high level of multiplicity in high-mass star formation, it is often difficult to disentangle different outflow contributions as described above for a few sources. We therefore used our outflow observations in context with other more detailed case studies for the sample to estimate the outflow directions. In almost all cases the velocity gradients seen in CH3CN (12$_3$–11$_3$) are roughly perpendicular to the direction of the molecular outflows, supporting the explanation that the velocity gradients probed by dense gas tracers are associated with rotational motions, therefore making these objects excellent disk candidates. The exceptions to this finding are W3IRS4 and IRAS 23033 MMS1a where the CH3CN velocity fields seem to trace the outflows, and G084 and G100 that have very weak velocity gradients. As described in the previous section, these four cores have been removed from our list as potential disk candidates. For the remainder of this work, we study the properties of the remaining 13 best disk candidates listed in Table 3. A plot showing the absolute difference between the assumed disk and outflow position angles is shown in Fig. C.7.
# 5. Analysis and discussion

## 5.1. Temperature distribution

We make use of the XCLASS software (Möller et al. 2017) within CASA to fit pixel-by-pixel the spectra of CH3CN (12K–11K) K = 0–6, including CH3¹³CN (12K–11K) K = 0–3, under the assumption of LTE, which is typically the case in such dense environments. In summary, the user provides initial guesses and sets allowed ranges for a set of parameters (column density, rotational temperature, line velocity, linewidth, and source size); the software solves the radiation transfer equation and through a minimisation routine modifies these parameters until the best fit is found. The software then creates maps for each of the fitted parameters. The details of the XCLASS fitting routine are outlined in Appendix B of Ahmadi et al. (2018) where we had fixed the source size to yield a beam filling factor of 1 and only fitted for the optically thin CH3CN (12K–11K) K = 4–6 lines, including CH3¹³CN (12K–11K) K = 0–3 isotopologues. While the observations presented here are resolved, our extensive tests with XCLASS have shown that having the source size parameter (and hence the filling factor) as a fitting parameter allows for much better simultaneous fits to the spectra when CH3CN (12K–11K) K = 0–3 lines are optically thick and the optically thinner CH3CN (12K–11K) K = 4–6 transitions<fnref n="5" />. It is important to note that the K = 3n transitions (i.e. 3 and 6) have double the statistical weight of the other transitions. For each source, we set the expected observed isotopic ratio according to its galactocentric distance, DGC, based on estimates of Wilson & Rood (1994): ¹²C/¹³C = (7.5 ± 1.9)DGC + (7.6 ± 12.9). Furthermore, we only fitted pixels above 6σ detection limit (σ values and isotopic ratios are listed in Table 3). The exceptional modelling approach for NGC 7538 IRS1 is summarised in Appendix D.1.

In Fig. 3, we show the rotational temperature maps. Assuming the lines are thermalised, rotational temperature can be directly taken as the kinetic temperature of the gas in these dense regions. Median temperature values are listed in Table 4 and are between 70–210 K with an uncertainty of 10–20% (see Appendix B of Gieser et al. 2019). It is useful to examine the maps of all fitted parameters, namely CH3CN column density, line velocity, linewidth, and source size, when analysing the spatial distribution of gas temperatures (see Appendix E). The XCLASS maps for the three cores that show CH3CN emission but that are not good disk candidates (G100, G084, W3IRS4) are shown in Fig. E.5. In the following, we summarise some of the trends we find for the physical conditions probed by CH3CN.

In the case of an idealised disk, the temperature and column density are expected to decrease with increasing radius. While for some of the sources in our sample we do see such behaviour (e.g. AFGL 2591, W3(H₂O)E, W3(H₂O)W, Cep A), the distribution of temperature is often not as smooth and symmetric as expected. Associating substructure in the temperature maps with underlying substructure in the disk is not straight forward. In fact, we have shown through simulated observations that the temperature is overestimated in regions of the disk with non-Gaussian line shapes, as well as multi-component line profiles to which we only fit one component (Ahmadi et al. 2019). Such multi-component line profiles are in fact expected as we start to resolve the infalling material from the envelope onto the disk, which results in gas motions at various velocities along the line of sight

<table><caption>Table 4. Gas mass estimates.</caption><thead><tr><th>Source</th><th>F<sub>v</sub><sup>(a)</sup><br>(mJy)</th><th>T<sub>med</sub><br>(K)</th><th>M<sub>gas</sub>(T<sub>med</sub>)<br>(M<sub>☉</sub>)</th><th>M<sub>gas</sub><br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>IRAS 23151</td><td>54</td><td>99</td><td>3.4</td><td>3.6</td></tr><tr><td>IRAS 23033 MMS 1b</td><td>47</td><td>84</td><td>3.9</td><td>3.9</td></tr><tr><td>AFGL2591</td><td>156</td><td>132</td><td>7.1</td><td>6.2</td></tr><tr><td>G75.78</td><td>123</td><td>102</td><td>9.7</td><td>8.1</td></tr><tr><td>IRAS 21078</td><td>301</td><td>120</td><td>3.2</td><td>3.3</td></tr><tr><td>G094</td><td>27</td><td>131</td><td>2.8</td><td>3.0</td></tr><tr><td>Cep A</td><td>457</td><td>212</td><td>0.8</td><td>0.6</td></tr><tr><td>NGC 7538 IRS9</td><td>55</td><td>212</td><td>0.9</td><td>1.2</td></tr><tr><td>W3(H<sub>2</sub>O)E</td><td>669</td><td>149</td><td>9.9</td><td>9.4</td></tr><tr><td>W3(H<sub>2</sub>O)W</td><td>544</td><td>148</td><td>8.1</td><td>7.0</td></tr><tr><td>IRAS 23385</td><td>96</td><td>72</td><td>18.3</td><td>19.1</td></tr><tr><td colspan="5"><strong>Pilot study</strong></td></tr><tr><td>NGC 7538 IRS1</td><td>2341<sup>(*)</sup></td><td>155</td><td>53.9</td><td>31.7</td></tr><tr><td>NGC 7538 S</td><td>67</td><td>132</td><td>2.0</td><td>1.8</td></tr></tbody></table>

Notes. (a) Continuum flux density at 1.37 mm within 6σ threshold (see text for details). (*) Subtracting the free-free contribution corresponds to a flux density of 2078 mJy at 1.37 mm (see text for details). Gas masses are calculated using the continuum emission maps with the free-free contribution subtracted.

depending on the inclination of the object. Furthermore, we do not account for density and temperature gradients along the line of sight.

For some of the sources we find the outskirts to be warmer than the inner regions (e.g. IRAS 23151, NGC 7538 IRS9). These regions have been properly modelled by XCLASS and are not the result of fitting noisy spectra as we have masked out the regions below 6σ (although the spectra are noisier in the outskirts than inside the cores). For these sources, it is helpful to simultaneously look at the maps of other parameters (see Figs. E.1–E.4). These outer hot components are often associated with large linewidths, consistent with our decomposed maps of linewidth (see Fig. B.3) that are able to cover more extended emission than the second moment maps. In many cases, the hot component associated with broad linewidths are found in the ejection direction of the outflows. This hints at the possibility that these regions are regions that may have been heated by the outflows and associated shocks and/or regions that have been carved out by the molecular outflows allowing for a deeper look into the cores. This touches on an important point about high-mass star formation, that such bipolar outflows that form early may be critical in creating cavities through which intense radiation pressure is able to escape (e.g. Yorke & Sonnhalter 2002; Krumholz et al. 2005; Kuiper et al. 2015, 2016; review by Pudritz & Ray 2019), allowing the accretion of mass onto their disks. The case of W3(H₂O) is perhaps the most impressive (see the large linewidths and positions of outflows in Fig. B.3).

## 5.2. Gas masses

Assuming optically thin dust emission at 1.3 mm and a gas-to-dust mass ratio R of 150 (Draine 2011), we can convert the flux density F<sub>v</sub> of the continuum observations to a mass via

$$M = \frac{d^2 F_v R}{B_v(T_D) \kappa_v}, \qquad (1)$$

<footnote n="5">⁵ See Fig. 5 of Mangum & Shirley (2015) to see the effect of beam filling factor on excitation temperatures.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Rotational temperature maps obtained by fitting CH<sub>3</sub>CN (12<sub>K</sub>–11<sub>K</sub>) *K* = 0–6 and CH<sub>3</sub><sup>13</sup>CN (12<sub>K</sub>–11<sub>K</sub>) *K* = 0–3 lines with XCLASS. The contours correspond to the 1.37 mm continuum as described in Fig. 1. For NGC 7538 IRS1, only the region outside the continuum to the south-west is modelled by XCLASS and was scaled towards the continuum peak position following a temperature power-law distribution *T* ∝ *r*<sup>-0.4</sup>. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel. Maps of column density, velocity offset, linewidth, and source size are shown in Appendix E.</figcaption></figure>
<table><caption>**Table 5.** *KeplerFit* parameters and dynamical mass estimates.</caption><thead><tr><th>Source</th><th><i>M</i><sub>PV</sub><br>(M<sub>☉</sub>)</th><th><i>v</i><sub>cont</sub><sup>(a)</sup><br>(km s<sup>-1</sup>)</th><th><i>v</i><sub>fitted</sub><sup>(b)</sup><br>(km s<sup>-1</sup>)</th><th><i>p</i><sub>fitted</sub><sup>(c)</sup><br>(au)</th><th><i>p</i><sub>masked</sub><sup>(d)</sup><br>(au)</th><th>Best fitted transition</th></tr></thead><tbody><tr><td>IRAS 23151</td><td>23 ± 2</td><td>-53.9</td><td>-54.8</td><td>120</td><td>-400 → 500</td><td>CH<sub>3</sub>CN (12<sub>3,4,5</sub>-11<sub>3,4,5</sub>) stacked</td></tr><tr><td>IRAS 23033 MMS 1b</td><td>14 ± 1</td><td>-54.4</td><td>-53.6</td><td>88</td><td>-400 → 400</td><td>CH<sub>3</sub>CN (12<sub>3</sub>-11<sub>3</sub>)</td></tr><tr><td>AFGL2591</td><td>23 ± 1</td><td>-5.6</td><td>-5.8</td><td>-623</td><td>-1000 → 1000</td><td>CH<sub>3</sub>CN (12<sub>3</sub>-11<sub>3</sub>)</td></tr><tr><td>G75.78</td><td>19 ± 3</td><td>0.7</td><td>0.0</td><td>-1078</td><td>-1800 → 2300</td><td>CH<sub>3</sub>CN (12<sub>5</sub>-11<sub>5</sub>)</td></tr><tr><td>IRAS 21078</td><td>13 ± 1</td><td>-6.6</td><td>-6.5</td><td>100</td><td>-250 → 250</td><td>CH<sub>3</sub>CN (12<sub>5</sub>-11<sub>5</sub>)</td></tr><tr><td>G094</td><td>20 ± 3</td><td>-42.9</td><td>-43.9</td><td>0</td><td>-500 → 1800</td><td>CH<sub>3</sub>CN (12<sub>3</sub>-11<sub>3</sub>)</td></tr><tr><td>Cep A</td><td>24 ± 1</td><td>-7.3</td><td>-8.3</td><td>10</td><td>-500 → 500</td><td>OCS (18-17)</td></tr><tr><td>NGC 7538 IRS9</td><td>25 ± 3</td><td>-56.7</td><td>-57.1</td><td>370</td><td>-200 → 700</td><td>HC<sub>3</sub>N (24-23)</td></tr><tr><td>W3(H<sub>2</sub>O) E<sup>(*)</sup></td><td>9 ± 1</td><td>-51.0</td><td>-51.3</td><td>0</td><td>-500 → 500</td><td>CH<sub>3</sub>CN (12<sub>4</sub>-11<sub>4</sub>)</td></tr><tr><td>W3(H<sub>2</sub>O) W<sup>(*)</sup></td><td>10 ± 1</td><td>-47.0</td><td>-48.0</td><td>99</td><td>-∞ → 200</td><td>CH<sub>3</sub>CN (12<sub>3</sub>-11<sub>3</sub>)</td></tr><tr><td>IRAS 23385</td><td>17 ± 2</td><td>-50.7</td><td>-51.0</td><td>1215</td><td>-1200 → 2000</td><td>HC<sub>3</sub>N (24-23)</td></tr><tr><td colspan="7"><strong>Pilot study</strong></td></tr><tr><td>NGC 7538 IRS1</td><td>22 ± 2</td><td>-56.0</td><td>-55.1</td><td>310</td><td>0 → 900</td><td>C<sub>2</sub>H<sub>5</sub>CN (25<sub>22,4</sub>-24<sub>22,3</sub>)</td></tr><tr><td>NGC 7538 S</td><td>12 ± 1</td><td>-60.8</td><td>-60.7</td><td>0</td><td>-∞ → 950</td><td>CH<sub>3</sub>CN (12<sub>5</sub>-11<sub>5</sub>)</td></tr></tbody></table>

**Notes.** (a) The velocity of the line at the position of the continuum peak. (b) The fitted central velocity in the PV diagram. The central velocity of the PV diagram was initiated at the velocity of the line at the position of the continuum peak ($p_{cont}$) and allowed to shift ±1 km s⁻¹ in the fitting routine (see dashed horizontal lines in Fig. 4). (c) The fitted central position in the PV diagram. The central position of the PV diagram was initiated at the continuum peak position and allowed to shift ± a synthesised beam in the fitting routine (see dashed vertical lines in Fig. 4). (d) Masked out (i.e. excluded) positions with respect to the fitted position, $p_{fitted}$ (see dotted vertical lines in Fig. 4). (e) Based on A-array observations presented in Ahmadi et al. (2018).

where *d* is the distance to the source and $κ_v = 0.9 \text{ cm}^2 \text{ g}^{-1}$ the dust absorption coefficient of thin ice mantles after $10^5$ yr of coagulation at a density of $10^6 \text{ cm}^{-3}$ (Ossenkopf & Henning 1994). Furthermore, $B_v(T_D)$ is the Planck function for blackbody radiation, which at the wavelength under investigation nearly follows the Rayleigh-Jeans law and is linearly dependent on the dust temperature $T_D$. We assume that gas and dust are in equilibrium and use the temperature maps presented in the previous section (Fig. 3) and the continuum maps (shown as contours in e.g. Fig. 1) to create maps of gas mass for each core. The only source with significant free-free contribution is NGC 7538 IRS1 (Beuther et al. 2017) that we subtract from our 1.37 mm continuum map (see Appendix D.2). Summing over all pixels with continuum emission above $6\sigma$ gives a final mass value for each core, listed in Table 4. We also list gas mass estimates calculated using the median temperature values, as it is common practice to use one temperature value when good thermometers are not observed. For this calculation, the flux density is calculated by converting the units of each map from Jy beam⁻¹ to Jy pixel⁻¹, masking out the continuum emission below $6\sigma$ as well as regions that do not have any detection for CH₃CN, and summing over the remaining pixels. For IRAS 23033 MMS1b the emission is integrated over a circular region encompassing this core only. For W3(H₂O) E and W, we simply define a vertical boundary between the two cores as presented in Ahmadi et al. (2018). It is reassuring that gas mass estimates calculated using median temperatures are, in almost all cases, very close to the estimates that use the full temperature maps. The only exception is NGC 7538 IRS1, which is the brightest continuum source in the sample. This large difference is due to the fact that gas temperatures are warm (> 200 K) where the bulk of continuum emission is found, while regions outside of the continuum emission are colder, lowering the value of the median temperature (155 K), which results in a higher-mass estimate.

## 5.3. Rotation curves and dynamical mass estimates

In Fig. 4 we show the position-velocity (PV) maps of best fitted transitions listed in Table 5 along cuts across the strongest velocity gradients as depicted by dotted lines in Fig. 2. The yellow curves correspond to Keplerian rotation profiles with $v(r) = \sqrt{GM/r}$ for enclosed masses listed in Table 5 (described further below). It is clear that in many of the sources the PV diagrams do not resemble differential rotational motion of the gas well, but rather mimic rigid-body-like rotation (i.e. $v(r) \propto r$). Furthermore, the existence of excess emission in the two quadrants opposite to where the rotational motions of the disk are expected hints to the presence of infalling material from the envelope (e.g. see models of Ohashi et al. 1997a and Tobin et al. 2012), outflowing gas, and/or coarse spatial resolution. While weak (or lack of) differential rotation could be an indicator that the CH₃CN emission may not be tracing disks in these observations, we have shown through synthetic observations in Ahmadi et al. (2019) that in fact PV diagrams of poorly resolved disks do often mimic rigid-body-like rotation. In the PV plots of poorly resolved disks, the molecular emission appears stretched in the position direction with excess infall contributions appearing in opposite quadrants as the disk and envelope material become blended. Therefore, using PV diagrams to rule out the existence of disks is a challenging task.

From the PV plots, we see that the rotation curves of some sources like AFGL2591, G75.78, IRAS 21078, and NGC 7538 IRS1 better resemble differential rotation of gas in a disk as some gas is found at high-velocities close to the central object at zero offset, with lower velocity gas further out. The PV plot of Cep A HW2 is particularly interesting, as it has a very extended nature, with a significant amount of emission present in the top-left quadrant associated with infall motions. In fact, using the 3D velocity field of 6.7 GHz CH₃OH masers, Sanna et al. (2017) show that the velocity field in Cep A HW2 is dominated
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Position–velocity (PV) plots for best fitted transitions listed in Table 5 along cuts in the direction of rotation as depicted by dotted lines in Fig. 2. The width of the cut is the size of a synthesised beam to increase the signal-to-noise ratio. The PV plots of W3(H₂O) E and W3(H₂O) W make use of the A-array only observations in order to gain more resolution and disentangle the two cores from the circumbinary material. PV plots for all other sources make use of the full dataset observed in the CORE survey. The contours correspond to the 6σ level increasing in steps of 24σ. Yellow lines show the Keplerian rotation curves for enclosed masses, *M*<sub>PV</sub>, listed in Table 5 for each source. The inner radii that were masked out (i.e. excluded) in the fitting routine are shown as vertical dotted lines. The fitted central velocities and positions are shown as dashed lines. The zero offset corresponds to the continuum peak position. The crosses in the bottom right corners correspond to the spatial and spectral resolutions.
by an infall component from a radius of 900 au (~1.3") down to
300 au (~0.4"), where a rotational component associated with
an accretion disk becomes dominant.

To determine protostellar masses, we fitted Keplerian curves to the 6σ edges of emission in the PV diagrams shown in Fig. 4 (the outermost contours) using the *KeplerFit*<fnref n="6" /> package (Bosco et al. 2019) in Python. This method was first introduced by Seifried et al. (2016) who showed that fitting the edge of emission in PV diagrams provides a more accurate estimate for the enclosed mass than other methods. We tested this method with simulated observations for our purposes in Ahmadi et al. (2019) where we showed that depending on how resolved disk-envelope systems are, fits to either the inner or outer parts of the PV diagrams may provide more accurate protostellar mass estimates. We further showed that fitting the PV plots of poorly resolved disk-envelope systems result in masses that are overestimated because emission is stretched over a larger area and high-velocity components of a disk that typically span small radii close to the central object are beam-smeared to lower velocities. For this reason, often velocities have a constant value at regions closest to the central object and have to be masked out in order to obtain better fits to the PV curves.

Using the *KeplerFit* code, we fitted Keplerian rotation profiles ($v(r) = \sqrt{GM/r}$) to the 6$\sigma$ edges of emission in the PV diagrams of CH$_3$CN (12$_K$−11$_K$) $K = 3–5$ transitions as they are not heavily blended with other (detected) lines, as well as HC$_3$N (24−23) and OCS (18−17). We also created and fitted PV plots of CH$_3$CN (12$_K$−11$_K$) $K = 3–5$ transitions stacked together, which is particularly useful for sources where the emission of individual CH$_3$CN transitions is weak. The central positions of the PV plots are set to the continuum peak positions and allowed to be fitted ± a synthesised beam in the fitting routine. The central velocities are set to the line velocities at the continuum peak positions and allowed to be fitted ±1 km s$^{-1}$ in the fitting routine. To obtain estimates on the fitted masses due to the choice of the edge of emission, we fix the best fitted central positions and central velocities and fit Keplerian rotation profiles to the 5$\sigma$ and 7$\sigma$ edges of emission in the PV diagrams. In Table 5 we present the fit parameters corresponding to the PV diagrams of the best fitted transitions chosen based on the goodness of the fit determined according to the lowest combined $\chi^2$ value and confirmed by eye. The reported masses are the error-weighted averages for the 5–7$\sigma$ solutions. The estimated dynamical masses, which for a Keplerian disk is the mass of the central object, are in the range ~10–25 $M_\odot$ for our sample of disk candidates. The Keplerian curves corresponding to the best mass estimates are shown as solid yellow lines on the PV diagrams shown in Fig. 4. The inner radii that were excluded in the fitting routine are shown as dotted vertical lines. The fitted central velocities and positions are also listed in Table 5 and are shown as dashed lines in Fig. 4.

The PV plot of NGC 7538 IRS1 shows the emission of CH3CN to be very extended in position, such that the best fitted PV curves yield masses that are extremely large and unrealistic considering the luminosity of the region. Observations of thermal CH3OH in the centimetre regime have shown that this source harbours two disk-like structures embedded within the same rotating circumbinary envelope (Beuther et al. 2017), consistent with earlier observations of multiple clusters of methanol masers at 6.7 GHz (Moscadelli & Goddi 2014). Therefore, fitting Keplerian curves to the PV diagrams of CH3CN, which likely

traces this circumbinary envelope, is not valid. For this reason,
the PV diagrams of the rarer molecule C2H5CN that better trace
regions closer to the protostars resulted in better fits. For the cir-
cumbinary envelope encompassing W3(H2O) E and W, we were
able to disentangle the gas contributions from each of the two
cores by using only the observations from the most extended
configuration (A-array) of NOEMA (Ahmadi et al. 2018). In
fact, the emission in the bottom-left quadrant of the PV plot for
W3(H2O) W corresponds to W3(H2O) E, therefore the expected
Keplerian curve corresponding to a mass of 10 M⊙, shown as
a yellow line, properly fits the contribution from W3(H2O) W
in the top-right quadrant. Similarly, NGC 7538 S likely harbours
multiple protostars as the continuum emission has an elongated
shape and the CH3CN (123-113) emission has multiple peaks
along this elongation (see Fig. A.1). Therefore, for this source,
only one side of the PV plot was fitted. Furthermore, the emis-
sion of CH3CN lines in IRAS 23385 lacks the high velocity
components expected in Keplerian rotation, and can be better
traced by HC3N(24-23) and OCS (18-17). In fact, thorough anal-
ysis of the CORE data by Cesaroni et al. (2019), revealed that the
velocity field of the main core in this source can be described by
a self-gravitating disk rotating about a ~9M⊙ star.

We compare the mass estimates from fitting the PV dia-
grams, *M*<sub>PV</sub>, against masses derived from the source luminosities
*M*<sub>lum</sub> using the spectral type-temperature calibration of Crowther
(2005, see Table 1). Masses derived from source luminosities
are in the range 19–41 *M*<sub>☉</sub>. Since the reported luminosities are
from multi-wavelength SED fitting of lower resolution observa-
tions (see Mottram et al. 2011), estimates for *M*<sub>lum</sub> are upper
limits, especially in sources with a high degree of fragmentation
and multiple cores detected (e.g. IRAS 23033 or W3(H<sub>2</sub>O)). Our
mass estimates from fitting the PV diagrams are either lower than
or consistent with these upper mass limits, within the estimated
errors. The reported dynamical masses and their uncertainties
should be used cautiously as they rely on the choice of fitted
molecular transition and whether the modelled transition traces
the disk well. We select the upper limit of the dynamical mass
estimates for the protostellar masses in the remainder of this
work. When the dynamical mass upper limits are higher than
the masses derived from the source luminosities, we select the
luminosity derived masses. For sources that were studied in more
depth using the same CORE data, we adopt the mass estimates
reported in these works (see Table 7 for a summary).

5.4. *Timescales*

Following Cesaroni et al. (2006), Beltrán et al. (2011) defined the term ‘toroids’ for rotating structures surrounding the most massive stars that span thousands of au with masses of several 10 $M_\odot$, comparable to or greater than the mass of the central protostar(s). These objects are thought to be transient, feeding a central cluster of forming stars, each surrounded by their own centrifugally supported disk on smaller scales (e.g. see Fig. 2 of Kuiper & Hosokawa 2018). They are analogous to rotating infalling envelopes around forming low-mass stars. Beltrán et al. (2011) studied the stability of rotating structures around O- and B-type stars by plotting the ratio of free-fall timescale to the rotational period at the outer radius, $t_{ff}/t_{rot}$, versus $M_{\text{gas}}$. The free-fall timescale represents the time it takes the gas to collapse under its own gravity. The rotational period at the outer radius is the time it takes the rotating structure to re-adjust its internal structure and stabilise after accreting new material. Simplifying the equations given in Sánchez-Monge et al. (2014), we calculate the rotational

<footnote n="6">⁶ Developed by Felix Bosco, available at https://github.com/felixbosco/KeplerFit</footnote>
period via

$$ 
\begin{aligned}
t_{\text{rot}} &= \frac{2\pi R}{v_{\text{rot}}} \\
&= 29.8 \left(\frac{R}{\text{au}}\right) \left(\frac{v_{\text{rot}}}{\text{km s}^{-1}}\right)^{-1} \text{yr},
\end{aligned} 
\quad (2) $$

where $v_{\text{rot}}$ is the rotational velocity at the outer radius, R. The free-fall timescale is calculated according to

$$ 
\begin{aligned}
t_{\text{ff}} &= \sqrt{\frac{3\pi}{32G\rho}} \\
&= 0.195 \left(\frac{M_{\text{gas}}}{M_{\odot}}\right)^{-0.5} \left(\frac{R}{\text{au}}\right)^{1.5} \text{yr}.
\end{aligned} 
\quad (3) $$

While this expression is valid if there is no central protostar (i.e. for a starless core), we use this expression to be consistent with the analysis presented in Beltrán & de Wit (2016).

Because the rotating structures traced by CH$_3$CN are not always symmetrically distributed on the position of the continuum peak, we calculate the rotational period and free-fall timescale for the redshifted and blueshifted sides of the rotating structures separately. We estimate the radii of the blueshifted and redshifted sides ($R_{\text{blue}}$ and $R_{\text{red}}$) as the distance from the continuum peak to the edge of the CH$_3$CN emission along the strongest velocity gradient (dotted lines on Fig. 2). We manually set the positions of $R_{\text{blue}}$ and $R_{\text{red}}$ if the velocity gradient appears to have a break close to the edge of the emitting region (e.g. blueshifted side of Cep A) or if the cores are blended (e.g. blueshifted side of W3(H$_2$O) W). The positions of $R_{\text{blue}}$ and $R_{\text{red}}$ are marked by an × in Fig. 2. $v_{\text{rot}}$ is estimated by subtracting the peak velocity at these positions from the LSR velocity measured at the peak of the continuum emission (listed in Table 3).

In Fig. 5, we show the plot of $t_{\text{ff}}/t_{\text{rot}}$ versus $M_{\text{gas}}$ for the redshifted and blueshifted sides of the 13 cores under investigation. The black dashed line corresponds to the theoretical curve of $t_{\text{ff}}/t_{\text{rot}}$ for a disk of mass $M_{\text{gas}}$ rotating about a star with mass $M_* = 10 M_{\odot}$. This is obtained by substituting $v_{\text{rot}} = \sqrt{G(M_{\text{gas}} + M_*)/R}$, that is the expected velocity at radius R if the rotation follows a Keplerian profile, in the above equations such that

$$ \frac{t_{\text{ff}}}{t_{\text{rot}}} = \sqrt{\frac{M_{\text{gas}} + M_*}{32 M_{\text{gas}}}}. \quad (4) $$

In a review by Beltrán & de Wit (2016), the authors created the $t_{\text{ff}}/t_{\text{rot}}$ versus $M_{\text{gas}}$ plot for the largest sample of intermediate and high-mass protostars from the literature and confirmed the findings of Beltrán et al. (2011) that disks and toroids are kinematically and dynamically different structures. Toroids contain more gas mass and are less rotationally supported with lower $t_{\text{ff}}/t_{\text{rot}}$, while for centrifugally supported disks in Keplerian rotation the infalling material has enough time to settle onto the disk, therefore the rotational timescale is shorter and the $t_{\text{ff}}/t_{\text{rot}}$ is higher. Comparing Fig. 5 to Fig. 14 of Beltrán & de Wit (2016), all of our sources lie in the regime corresponding to centrifugally supported disks, while toroids reside in a regime below the theoretical curve where gas masses are larger than 40 $M_{\odot}$, beyond the limits plotted in Fig. 5.

## 5.5. Specific angular momentum

For a molecular cloud to form a star, it must lose many orders of magnitude in specific angular momentum (j). Various mechanisms at a multitude of scales can contribute to this loss

<figure><img src="image_6.png" /><figcaption>Fig. 5. Ratio of free-fall to rotational timescale as a function of disk gas mass, showing most of the sources with CH$_3$CN velocity gradients are rotationally supported. The blue and red dots correspond to the values calculated based on the rotational velocities at the edges of the disks on the blueshifted and redshifted sides, as depicted by crosses in Fig. 2 and listed in Table 3. The black dashed line corresponds to theoretical curve of $t_{\text{ff}}/t_{\text{rot}}$ for a disk of mass $M_{\text{gas}}$ in Keplerian rotation about a star with mass $M_* = 10 M_{\odot}$ (see Eq. (4)).</figcaption></figure>

(and transport) of angular momentum (e.g. disks, jets, outflows, winds, magnetic fields). Specific angular momenta have predominantly been measured in regions of low-mass star formation. In particular, it has been shown that on dense core scales (0.1 pc), the local specific angular momentum (i.e. the product of $v_{\text{rot}}$ and r) scales as $r^{1.6}$ (Goodman et al. 1993), down to a scale of about ~6000 au where the specific angular momentum is constant down to disk scales (~100 au). Ohashi et al. (1997b) was first to measure a constant j value of ~$10^{-3}$ km s$^{-1}$ pc for this intermediate regime and attribute it to protostellar envelopes that are dynamically collapsing. Below ~100 au, Keplerian rotation dominates and once the disk becomes unstable, material begins to move inwards and angular momentum outwards. A review by Belloche (2013) compiles most of the relevant observations from literature on this topic and presents the correlation of the specific angular momentum with radius from disk to core scales in their Fig. 8. In a recent work, Pineda et al. (2019) calculated the specific angular momentum radial profiles towards two Class 0 objects and a hydrostatic core candidate in Perseus and found a power-law relation with $j \propto r^{1.8}$ from 800–10 000 au. The authors concluded that because the power-law index is between that of solid-body rotation ($j \propto r^2$) and pure turbulence ($j \propto r^{1.5}$), the influence of the dense core’s initial level of turbulence is still present even down to 1000 au, inherited from the turbulent parental molecular cloud. Furthermore, Gaudel et al. (2020) calculated the specific angular momentum profiles of 11 Class 0 protostellar envelopes from the CALYPSO survey from ~50–5000 au scales and found a decreasing power-law relation with $j \propto r^{1.6 \pm 0.2}$ down to ~1600 au and relatively constant down to ~50 au. Since these findings are based on observational evidence from low-mass star-forming regions, it is important to calculate the specific angular momenta for our more massive sources for comparison. Additionally, it is important to note that the same scaling relations may not apply in high-mass star-forming regions as it remains unclear whether high-mass stars form via ordered collapse of pre-existing cores or in a more dynamic environment.
In an effort to determine whether the rotational signatures that we have observed are more in agreement with core rotation, conservation of angular momentum, or differential rotation of a disk-like object, we calculate the specific angular momenta as a function of radius following the prescription in Belloche (2013). The rotational kinetic energy of an object in solid body rotation with mass M and radius r is $E_{rot} = \frac{1}{2}I\Omega^2$, where $\Omega$ is the angular velocity and $I = \frac{2}{3}Mr^2(\frac{3-\alpha}{5-\alpha})$ is the moment of inertia (Belloche 2013). The value of $\alpha$ is determined based on the power law density profile of the object such that $\rho \propto r^{-\alpha}$. Gieser et al. (2021) fitted the continuum visibility profiles of the sources in the CORE sample and found a mean density power law index of $2 \pm 0.2$. Therefore, setting $\alpha = 2$ and knowing the angular momentum $J = I\Omega$, the specific angular momentum is

$$j = \frac{J}{M} = \frac{2}{9}\Omega r^2. \qquad (5)$$

The angular velocity is calculated by subtracting the LSR velocity, defined to be the velocity at the position of the continuum peak (see Table 3), from the centroid velocity $v_{cen}$, at each position and dividing by the radius at that position with a $\sin(i)$ correction factor included for the unknown inclination $i$:

$$\Omega = \frac{|v_{cen} - v_{LSR}|}{r \sin(i)}. \qquad (6)$$

Putting together Eqs. (5) and (6), we can calculate the specific angular momentum as a function of radius times the unknown $\sin(i)$ correction for each of the sources in the sample, such that

$$j_{obs}(r) \sin(i) = \frac{2}{9} |v_{cen} - v_{LSR}| r. \qquad (7)$$

The observed specific angular momentum for the CORE sample were calculated using the decomposed velocity maps of CH$_3$CN (12$_3$–11$_3$) for each source (see Fig. B.2) along the cut with the strongest velocity gradient (dotted lines) for the redshifted and blueshifted sides separately. The resulting radial profiles are presented in Fig. 6 as circular data points with their errors, propagated from the error in velocity from the Gaussian decomposition, shown as regions with filled colours. The upward and downward pointing triangles correspond to the same procedure but performed for a cut with a position angle that is respectively 10° larger and smaller than the position angle of the cut with the strongest velocity gradient. The distance between two data points corresponds to half-beam spacing to ensure the data points are independent of each other, and the specific angular momenta plotted are averaged over a region of 3 × 3 pixels. Only the sources for which we had at least four data points (i.e. more than two beams across) for either the redshifted or blueshifted sides are shown. Therefore, this criterion excludes the core IRAS 23385.

We can compare the angular momentum of each source with some estimate for the maximum value above which the core becomes gravitationally unbound. Because neither the rotation nor the density profile is expected to be uniform, we follow the approach of Zhang et al. (2018), assuming the structure is flattened with its mass growing linearly with radius (see Mestel disk: Mestel 1963). The maximum specific angular momentum can be written as

$$j_{max}(r) = \frac{1}{2} \sqrt{G M(r) r}, \qquad (8)$$

where $M(r)$ is the mass interior to $r$ and is equal to $M_* + M_{gas}(r)$. For the protostellar masses, we used the estimates from fitting the PV diagrams of dense gas tracers (see Sect. 5.3) unless other detailed studies of the sources using our data provided better estimates (see Table 7). We used the gas mass maps created using Eq. (1) and the temperature maps obtained from our radiative transfer modelling with XCLASS to calculate $M(r)$. The maximum specific angular momenta are over-plotted on Fig. 6 with black lines, and the grey shaded areas surrounding the black curves correspond to the estimate for $j_{max}$ with $M_* \pm 5 M_{\odot}$. The maximum radius for which we calculate $j$ shown in blue and red data points corresponds to the extent of the blueshifted and redshifted gas emission, while the extent of the black curves corresponding to $j_{max}$ mark the edge of the continuum emission at $6\sigma$. Comparing this curve with the observed $j$ radial profiles, we see that almost all cores are gravitationally bound as they lie below $j_{max}$ with the exception of NGC 7538 S which has an elongated shape in both the continuum and CH$_3$CN towards the redshifted side – a sign that further fragmentation has taken place.

From the specific angular momentum radial profiles (see Fig. 6) we see that all sources span a range in $j$ of 1–2 orders of magnitude, on average around ~$10^{-3}$ km s$^{-1}$ pc, and exhibit a shallowing of their power-law slopes with increasing radius. The dashed lines correspond to least-squared power-law fits to the first two datapoints and extrapolated to larger radii to clearly show this behaviour. The power-law slopes and intercepts corresponding to these dashed lines are listed in Table 6, with a median slope of $j \propto r^{1.7}$. For some sources (e.g. IRAS 23151, IRAS 23033 MMS1b, G094, NGC 7538 S), the flattening of $j$ on larger scales is more subtle than other sources (e.g. G75.78 and Cep A). For AFGL2591 and G75.78 the decrease in $j$ is followed by a sharp increase as a result of the complex velocity structures in the outskirt of these objects (see Fig. 2), likely as a result of material flowing from larger scales. The redshifted side of W3(H$_2$O) E and the blueshifted side of W3(H$_2$O) W are contaminated by the emission from each other and have a similar shape, while the blueshifted side of W3(H$_2$O) E and the redshifted side of W3(H$_2$O) W which actually correspond to the edges of these rotating structures show the decrease in $j$ more clearly.

In order to better understand the observational findings, we analyse the expected specific angular momentum profiles for the numerical simulations presented in Oliva & Kuiper (2020) as well as synthetic observations of the system at varying inclinations and angular resolutions to study their effects on $j$ (see Appendix F). We have purposefully not corrected for the inclinations of the simulated observations in order to compare the findings with our CORE observations with unknown inclinations. Similar to the observational findings for low-mass star formation summarised in the first paragraph of this section, we find three distinct regimes in the $j$ radial profile of the simulations (see Fig. F.1). The inner disk region ($r < 100$ au) actively accreting material onto the central protostar with $j \propto r^{0.6}$ close to Keplerian rotation ($j \propto r^{0.5}$), an intermediate region where $j$ is quite flat but not exactly constant, and the region beyond ~500 au is rather consistent with the $j \propto r^{1.6}$ relation that had been found for the rotation of low-mass dense cores and Class 0 envelopes (e.g. Gaudel et al. 2020). As resolution worsens, the substructure in the $j$ radial profile becomes smoother with the power-law relations for the different regimes becoming more continuous and not as clear to disentangle. At the scale of synthetic NOEMA observations at 2 kpc, the envelope and disk components become completely blended, with only one power-law relation seen across all scales, with $j \propto r^{1.6}$. This is
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Specific angular momentum radial profiles calculated using Eq. (7) along the cut with the strongest velocity gradient (dotted lines in Fig. 4), (circles), a cut with position angle +10° (triangles pointing up), and a cut with position angle -10° (triangles pointing down) with respect to the strongest velocity gradient, for the redshifted (red) and blueshifted (blue) sides of the emission. The blue and red filled regions show the error in $j \sin(i)$ for the circles. We note that the inclinations are not corrected for as they are not known. The blue and red dashed lines show least-squared power-law fits to the first two data points and extrapolated to larger radii. The values of the fitted slopes are noted in the legends and listed in Table 6 along with the intercepts. The blue and red dotted vertical lines correspond to radii associated with the edges of the disks depicted with × markers in Fig. 2. The black curve shows the maximum specific angular momentum calculated using Eq. (8) for $M_*$ listed in the sub-panels with the grey region corresponding to $M_* \pm 5 M_\odot$. The grey dashed line shows the $j \sin(i) \propto r^{1.7}$ slope.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. continued.</figcaption></figure>

consistent with our findings for the CORE sample with a median $j \propto r^{1.7}$ (see Fig. 6). It is also evident that at the spatial resolution of our observations, the lack of information on the inclination is not an important factor and does not impact our findings.

There exists a small shallowing of the slope (a kink) in the *j* profiles of poorly resolved simulated observations and the CORE observations. If the decrease in *j* is in fact due to the transition

between a Keplerian disk and the collapse of the envelope for which the angular momentum is constant (i.e. *j* is constant), one could put a limit on the size of the inner disk that is actively accreting onto the protostar. We refrain from calculating exact numbers for each source due to the uncertainties in such an estimate and our limited resolution, but note that for most of the sources this transition is seen around 1000–2000 au. The
<table><caption>**Table 6.** Fit parameters to the specific angular momentum radial profiles shown in Fig. 6 (dashed lines).</caption><thead><tr><th>Source</th><th>a<sub>blue</sub></th><th>a<sub>red</sub></th><th>b<sub>blue</sub></th><th>b<sub>red</sub></th></tr></thead><tbody><tr><td>IRAS 23151</td><td>1.8</td><td>1.4</td><td>-2.4</td><td>-3.0</td></tr><tr><td>IRAS 23033 MMS 1b</td><td>(...)</td><td>2.1</td><td>(...)</td><td>-2.8</td></tr><tr><td>AFGL2591</td><td>1.4</td><td>1.9</td><td>-3.3</td><td>-2.9</td></tr><tr><td>G75.78</td><td>1.5</td><td>1.6</td><td>-2.7</td><td>-2.7</td></tr><tr><td>IRAS 21078</td><td>1.3</td><td>1.6</td><td>-2.8</td><td>-2.5</td></tr><tr><td>G094</td><td>1.8</td><td>1.0</td><td>-3.0</td><td>-3.2</td></tr><tr><td>Cep A</td><td>1.8</td><td>1.6</td><td>-2.3</td><td>-2.5</td></tr><tr><td>NGC 7538 IRS9</td><td>1.3</td><td>1.5</td><td>-2.5</td><td>-2.4</td></tr><tr><td>W3(H<sub>2</sub>O) E</td><td>2.1</td><td>2.1</td><td>-3.2</td><td>-2.8</td></tr><tr><td>W3(H<sub>2</sub>O) W</td><td>1.7</td><td>1.9</td><td>-2.7</td><td>-2.5</td></tr><tr><td colspan="5"><hr></td></tr><tr><td>Pilot study</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 7538 IRS1</td><td>1.6</td><td>2.2</td><td>-2.3</td><td>-2.6</td></tr><tr><td>NGC 7538 S</td><td>3.7</td><td>1.7</td><td>-2.0</td><td>-2.3</td></tr></tbody></table>

**Notes.** The heading parameters *a* and *b* satisfy the following relation: $j \sin(i) = 10^b (\frac{R}{1000 \text{ au}})^a \text{ km s}^{-1} \text{ pc}$, for the blueshifted and redshifted sides.

only exceptions are Cep A and IRAS 21078 that are the closest sources in the sample and hence resolved the best, for which this transition is seen around 300–400 au, and 500–700 au, respectively.

## 5.6. Toomre stability

We determine the stability of the disk candidates against axisymmetric gravitational collapse using the Toomre *Q* parameter (Toomre 1964),

$$Q = \frac{c_s \Omega}{\pi G \Sigma} \qquad (9)$$

The Toomre equation determines the balance between the stabilising force of thermal pressure (through the sound speed, $c_s$) and shear forces (through the angular velocity, $\Omega$) against the local gravity of the disk (through the disk surface density, $\Sigma$). For a disk that is sufficiently massive, the self-gravity of the disk could cause it to fragment. The Toomre *Q* parameter quantifies this scenario such that if $Q < 1$, the disk is prone to fragmentation. While this critical threshold has been derived for axisymmetric instabilities in an infinitely thin disk, non-axisymmetric instabilities can form in thicker disks at higher values of *Q* (discussed further below and in Sect. 5.6.4).

The sound speed is calculated using the temperature maps obtained from our radiative transfer modelling of CH$_3$CN lines with XCLASS, presented in Sect. 5.1, such that

$$c_s = \sqrt{\frac{\gamma k_B T}{\mu m_H}}, \qquad (10)$$

where $\gamma$ is the adiabatic index with a value of 7/5, $k_B$ is the Boltzmann constant, $\mu$ is the mean molecular weight with a value of 2.8 (Kauffmann et al. 2008), and $m_H$ is the mass of the hydrogen atom. Using the obtained temperature maps and assuming gas and dust temperatures are coupled, the surface density of the disk is calculated according to

$$\Sigma = \frac{S_v R}{B_v(T_D) \theta_B \kappa_v}, \qquad (11)$$

where $S_v$ is the continuum intensity and $\theta_B$ is the beam solid angle. Furthermore, assuming the disk is in gravito-centrifugal equilibrium following a Keplerian rotation profile, the angular velocity at a given radius is

$$\Omega(r) = \sqrt{\frac{G(M_* + M_{\text{disk}}(r))}{r^3}}, \qquad (12)$$

where $M_*$ is the protostellar mass and $M_{\text{disk}}(r)$ is the gas mass interior to $r$ calculated from the continuum maps (see Sect. 5.2 and Table 7). Protostellar masses were estimated from fitting the PV diagrams of CH$_3$CN transitions and other dense gas tracers as discussed in Sect. 5.3. We select the upper limit of the dynamical mass estimates for the protostellar masses to calculate upper limits for the Toomre *Q* parameter. When the dynamical mass upper limits are higher than the masses derived from the source luminosities, we select the luminosity derived masses. For sources that were studied in more depth using the same CORE data (IRAS 23033: Bosco et al. 2019; IRAS 23385: Cesaroni et al. 2019), we adopt the mass estimates reported by these authors. All mass estimates and the chosen protostellar masses are summarised in Table 7.

Assuming the protostar is positioned at the location of the continuum peak, we calculate the Toomre *Q* parameter for each disk candidate, accounting for the self-gravity of the disk by including the enclosed gas mass at each radius. The resulting *Q* maps are presented in Fig. 7 for the 13 disk candidates in the CORE survey. Regions outside of the $6\sigma$ continuum contours are masked out since the continuum maps are used as input for the calculation of *Q*. The minimum and median *Q* values for each source are provided in Table 7.

The critical value of the Toomre parameter, $Q_{\text{crit}}$, which determines the stability of a disk is in theory equal to 1, while effectively it lies somewhere between 1–2 (Kennicutt 1989; Elmegreen 1991, 1994). Considering the uncertainties in our analysis, especially in determining temperatures from radiative transfer modelling of CH$_3$CN K-ladder (10–20% uncertainty, see Appendix B of Gieser et al. 2019) as well as the fact that the derived temperatures may be probing layers above the disk mid-plane and warmer regions shocked by outflows, we assume that regions where $Q \lesssim 2$ are unstable against gravitational collapse. Furthermore, it is important to keep in mind that this Toomre analysis depends on the spatial resolution and gives average values on size-scales of the beam. Toomre-stable disks could have smaller-scale substructures that may be gravitationally unstable which only higher spatial resolution observations can disentangle. We studied these caveats in depth by creating synthetic NOEMA observations of a numerical simulation of a gravitationally unstable and fragmented disk and confirmed the $Q_{\text{crit}} \approx 2$ to be a reasonable threshold for our CORE observations presented here (Ahmadi et al. 2019).

Because the Toomre parameter quantifies local stability in a disk, the global median values listed in Table 7 should be used with caution as a determining factor for whether a disk is stable against gravitational instabilities. In fact, numerical simulation of a massive fragmenting disk has shown that the disk can be globally Toomre-stable while locally Toomre-unstable and undergoing gravitational collapse (Oliva & Kuiper 2020). Hence having observations with sufficient spatial resolution and analysing the spatial distribution of the Toomre *Q* parameter rather than calculating a global value are crucial in studying the stability of such disks. In order to understand the local variations of *Q*, Fig. 8 shows the distribution of median and minimum *Q* values as a function of radius. While the global
<table><caption>**Table 7.** Overview of mass estimates and Toomre *Q* results.</caption><thead><tr><th>Source</th><th><i>M</i><sub>gas</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>PV</sub><sup>(a)</sup><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>lum</sub><sup>(b)</sup><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>*</sub><sup>(c)</sup><br>(<i>M</i><sub>☉</sub>)</th><th><i>Q</i><sub>min</sub></th><th><i>Q</i><sub>median</sub></th><th>(<i>H</i>/<i>r</i>)<sup>(d)</sup><sub>median</sub></th><th>Toomre stable</th></tr></thead><tbody><tr><td>IRAS 23151</td><td>3.6</td><td>23 ± 2</td><td>24</td><td>24</td><td>1.2</td><td>4.4</td><td>0.16</td><td>No</td></tr><tr><td>IRAS 23033 MMS 1b</td><td>3.9</td><td>14 ± 1</td><td>21</td><td>19<sup>(*)</sup></td><td>0.5</td><td>2.0</td><td>0.15</td><td>No</td></tr><tr><td>AFGL2591</td><td>6.2</td><td>23 ± 1</td><td>40</td><td>24</td><td>1.2</td><td>2.7</td><td>0.21</td><td>No</td></tr><tr><td>G75.78</td><td>8.1</td><td>19 ± 3</td><td>34</td><td>22</td><td>0.3</td><td>2.1</td><td>0.23</td><td>No</td></tr><tr><td>IRAS 21078</td><td>3.3</td><td>13 ± 1</td><td>19</td><td>14</td><td>0.4</td><td>3.3</td><td>0.23</td><td>No</td></tr><tr><td>G094</td><td>3.0</td><td>20 ± 3</td><td>25</td><td>23</td><td>0.8</td><td>6.4</td><td>0.23</td><td>No</td></tr><tr><td>Cep A</td><td>0.6</td><td>24 ± 1</td><td>21</td><td>21</td><td>5.1</td><td>11.9</td><td>0.17</td><td>Yes</td></tr><tr><td>NGC 7538 IRS9</td><td>1.2</td><td>25 ± 3</td><td>24</td><td>24</td><td>4.4</td><td>15.0</td><td>0.24</td><td>Yes</td></tr><tr><td>W3(H<sub>2</sub>O) E</td><td>9.4</td><td>9 ± 1</td><td>33</td><td>10<sup>(†)</sup></td><td>0.1</td><td>0.9</td><td>0.23</td><td>No</td></tr><tr><td>W3(H<sub>2</sub>O) W</td><td>7.0</td><td>10 ± 1</td><td>33</td><td>11<sup>(†)</sup></td><td>0.1</td><td>1.0</td><td>0.23</td><td>No</td></tr><tr><td>IRAS 23385</td><td>19.1</td><td>17 ± 2</td><td>21</td><td>9<sup>(‡)</sup></td><td>0.1</td><td>0.9</td><td>0.25</td><td>No</td></tr><tr><td colspan="9"><strong>Pilot study</strong></td></tr><tr><td>NGC 7538 IRS1</td><td>31.7</td><td>22 ± 2</td><td>41</td><td>24</td><td>0.3</td><td>0.7</td><td>0.16</td><td>No</td></tr><tr><td>NGC 7538 S</td><td>1.8</td><td>12 ± 1</td><td>21</td><td>13</td><td>0.9</td><td>4.2</td><td>0.26</td><td>No</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Dynamical mass estimates from fitting Keplerian curves to the $6\sigma$ edges of emission in the PV diagrams shown in Fig. 4 (see Sect. 5.3.) <sup>(b)</sup>Luminosity mass extrapolated from Table 1 of Crowther (2005) using luminosities reported in Table 1. <sup>(c)</sup>Protostellar mass used in the Toomre *Q* analysis. <sup>(d)</sup>Median ratio of disk scale height to radius. <sup>(e)</sup>Based on Bosco et al. (2019) using data from the CORE survey. <sup>(f)</sup>Based on A-array observations presented in Ahmadi et al. (2018) using data from the CORE survey. <sup>(g)</sup>Based on Cesaroni et al. (2019) using data from the CORE survey.

median *Q* value lies above the critical *Q* threshold for a number of sources, the minimum *Q* dips below the threshold in the outskirts of these disk candidates for all sources but two (Cep A and NGC 7538 IRS9). Therefore, we find 11 out of 13 disks to be Toomre unstable against gravitational collapse.

### 5.6.1. Toomre unstable disks

There exist regions with $Q \lesssim 2$ for most disk candidates (see Fig. 7), meaning they are unstable against gravitational instabilities and prone to fragmentation. While some candidate disks appear to be Toomre unstable almost entirely (e.g. G75.78, W3(H<sub>2</sub>O) E and W, IRAS 23385, NGC 7838 IRS1), others show signs of instability only in their outskirts (e.g. NGC 7538 S). For some sources, the low *Q* values are coincident with the elon- gation of the continuum, hinting at the existence of underlying substructure. The best examples for such cases are IRAS 23385, and NGC 7538 S. In fact, the in-depth analysis of Cesaroni et al. (2019) using the CORE survey data confirms that IRAS 23385 is forming a cluster of massive stars (marked by green plus symbols in Fig. 7). Another interesting case worth noting is IRAS 21078, for which the low *Q* values are coincident with a secondary source detected in continuum, ~0.5" to the south-east of the continuum peak.

Similarly, other sources for which we know further frag- mentation has taken place are AFGL2591 and NGC 7538 IRS1, marked by green plus symbols in Fig. 7. For AFGL2591, NOEMA observations at 843 µm (0.19" ~ 700 au – roughly a factor of 2 better than the work presented here and providing the highest achievable spatial resolution at (sub-)mm wavelengths in the northern sky) revealed multiple fragments within the inner ~1000 au, supporting the disk-fragmentation scenario (Suri et al. 2021). For NGC 7538 IRS1, Beuther et al. (2017) observed the fragmentation in centimetre continuum emission (with 0.055" resolution, ~150 au) and found two HCH II regions with a pro- jected separation of ~430 au. The authors resolved velocity

gradients in thermal CH<sub>3</sub>OH across each HCH II region, painting a picture of two disk-like structures embedded within the same rotating circumbinary envelope. Methanol masers at 6.7 GHz have also been detected, consistent with the accretion disk sce- nario (Moscadelli & Goddi 2014). The position angles of these disks are misaligned from that of the rotating structure we detect by about 30°. This scenario is similar to the synthetic obser- vations of the numerical simulation we presented in Ahmadi et al. (2019) whereby the disk on large scales has a different rotation axis than the inner small-scale disks that have formed around each of the fragments. Whether the large-scale rotating structure is a circumbinary envelope rather than a flattened disk- like structure is difficult to establish, especially because our line observations at the location of the protostar(s) are seen in absorp- tion. The large extent of CH<sub>3</sub>CN emission in the PV diagrams (see Fig. 4) at first glance may hint at the possibility that the structure is not a disk, especially because the mass estimates from fitting the PV diagram are much in excess of the mass esti- mated from the luminosity. However, as we showed in Ahmadi et al. (2019), PV diagrams of poorly resolved disks often have very large sizes resulting in mass estimates that are significantly higher than the true enclosed mass. Due to the small separation between the two embedded disks, it is likely that they formed via disk fragmentation.

### 5.6.2. Toomre stable disks

The candidate disks that appear to be completely Toomre sta- ble against gravitational instabilities include Cep A HW2 and NGC 7538 IRS9. The stability of Cep A is in part a result of warmer temperatures across the entire structure as compared with other sources as well as the low gas mass (0.5 $M_☉$) to stel- lar mass (20 $M_☉$) ratio. Therefore, the self-gravity of the disk is not very strong and the disk is warm prohibiting any possi- ble fragments to effectively cool and survive. The *Q* values for Cep A are so high that even if we have overestimated the mass
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Toomre *Q* maps for the best disk candidates in the CORE survey assuming a protostar is located at the position of the continuum peak as depicted by a star and accounting for the self-gravity of the disk (protostellar and gas mass values are listed in Table 7). The contours correspond to the 1.37 mm continuum as described in Fig. 1. Regions outside of $6\sigma$ continuum contours are masked out. Fragments detected by ancillary observations are marked by plus symbols in green for AFGL2591 (Suri et al. 2021), IRAS 23385 (Cesaroni et al. 2019), and NGC 7538 IRS1 (Beuther et al. 2017). The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Distribution of median and minimum *Q* as a function of radius shown in solid and dashed blue lines. The dotted horizontal line corresponds to the global median *Q* computed over the entire disk and listed in Table 7. The dash-dotted red horizontal line shows the critical *Q* = 2 threshold.</figcaption></figure>

of the protostar by a factor of 2 (e.g. due to contamination of the blueshifted emission by a hot core that is located 0.3" away from the main continuum peak (Jiménez-Serra et al. 2009), the disk would still appear stable at the scales probed by our observations. The low gas to stellar mass ratio of these two systems could be an indication that these systems may be somewhat closer to the end of their accretion phase, though Cep A HW2 is considered IR-dark while NGC 7538 IRS9 is IR-bright.

### 5.6.3. General trends

In Fig. 9 we present the *Q* distribution versus the ratio of gas to stellar mass, coloured according to the YSO's luminosity. The boxes extend from the first to the third quartiles, with a red line at the median *Q* value (also listed in Table 7). As explained in the previous section, since Toomre stability is a measure for the *local* stability of the disk against axisymmetric gravitational instabilities (see also e.g. Oliva & Kuiper 2020), simply using $Q_{crit} \approx 2$ to determine the stability of these disks from the median *Q* values is incorrect. Based on the analysis presented above and by looking at the distribution of *Q* for each individual case against $Q_{crit}$ (see Figs. 7–9), disks with masses greater than ~10–20% of the mass of their host (proto)stars are subject to fragmentation through gravitational instabilities. Numerical experiments of Kratter et al. (2010) find that the maximum disk mass before disk fragmentation occurs is on the same order as the stellar mass, with a review on the subject by Kratter & Lodato (2016) providing a value of $M_{disk}/M_* \ge 0.1$ for gravitational instabilities to become important. Our observations confirm that disk fragmentation must be prevalent in high-mass star formation as

we find 11 of 13 disks studied here to be gravitationally unstable, all of which have $M_{gas}/M_* \ge 0.1$. Furthermore, we find an inverse relationship between *Q* and $M_{gas}/M_*$, a result of the expected tight anti-correlation between *Q* and $M_{gas}$ and less so due to a correlation between *Q* and $M_*$ (see Fig. 10). We also find indications that more luminous YSOs tend to have disks that are more massive and as a result more prone to fragmentation (see left panel of Fig. 10). This could be a selection effect as high-mass protostars with lower-mass disks must exist at a later stage of accretion. Additionally, the luminosities are derived from lower resolution observations corresponding to a different number of cores in each star-forming region.

Oliva & Kuiper (2020) investigated the formation and early evolution of an accretion disk surrounding a forming high-mass protostar in their 3D radiation-hydrodynamic simulation with high enough resolution to resolve the pressure scale height and Jeans length of the disk. The accretion disk with Keplerian-like rotation is globally Toomre-stable but locally Toomre-unstable and therefore fragments. The fragments are dynamic within the system, migrating inwards and can even get flung out of the system if they are on highly eccentric orbits. More massive fragments that remain on more stable orbits form their own accretion disks and become hydrostatic cores, some of which go on to form second Larson cores. The authors further study the fate of these companions and conclude that the hydrostatic core fragments migrate inwards and accrete onto the central protostar causing accretion bursts, while the second cores would likely become spectroscopic companions to the protostar. They find that companions produced through disk fragmentation can exist at a range of distances from the primary protostar, from spectroscopic
<figure><img src="image_11.png" /><figcaption>Fig. 9. Box plot showing the Q distribution versus the ratio of gas to stellar mass, coloured according to the luminosity of the regions within which they reside. The boxes extend from the first to the third quartiles, with a red line at the median Q value. The whiskers extend from the box by ±1.5 times the box size. Disks with Q ≲ 2 are Toomre unstable (i.e. all but Cep A & NGC 7538 IRS9).</figcaption></figure>

multiples to companions at distances between 1000 and 2000 au. Synthetic ALMA and NOEMA observations of this numerical simulation are presented in Ahmadi et al. (2019).

Whether fragments survive as companions to the central protostar depends largely on whether they can cool fast enough. The cooling time has been confirmed to be short enough relative to the orbital period both numerically (Klassen et al. 2016; Meyer et al. 2018) and observationally (Ahmadi et al. 2018) such that fragments that are beyond ~50 au of the central protostar should be rapidly cooling (see review by Kratter & Lodato 2016). This is an important finding as nearly all high-mass stars are found in binary or multiple systems (Chini et al. 2012; Sana et al. 2012, 2014). With evidence that massive binaries start their lives in wide pairs and harden into tighter massive binaries over time (Ramírez-Tannus et al. 2021), disk fragmentation may be a possible mechanism creating the initial binary and multiple systems at large separations before the inward migration occurs. Additionally, disk fragmentation can have profound implications for the formation of high-mass stars, as the mass flow is shared by multiple accretors rather than a single protostar (e.g. Bate & Bonnell 1997; Bate 2000; Peters et al. 2010b).

It is important to note that we have neglected the effect of magnetic fields as we do not have information about the magnetic field strengths in our survey. Theoretically, one can incorporate the Alfvén velocity in the Toomre equation in order to determine the magnetic Q parameter, Qmag (Kim & Ostriker 2001). In their numerical collapse simulations of strongly magnetised and turbulent molecular cloud cores, Seifried et al. (2013) compare Q measurements with Qmag. They find that Qmag rises above the value of Q, confirming that disk fragmentation is suppressed by the presence of magnetic fields as magnetic pressure provides added support against collapse. However, fragmentation in a magnetised environment can still proceed through other mechanisms (e.g. high amplitude fluctuations in the core initially, ambipolar diffusion and fragmentation during the second collapse; see e.g. Hennebelle & Teyssier 2008). In fact, if the disk is well within the non-ideal MHD limit, the circumstellar disk structure is very similar to the non-magnetic case (e.g.

Kölligan & Kuiper 2018). Numerical simulations of massive protostellar disks that include non-ideal MHD and radiative transfer have shown that magnetic pressure is not as supportive as thermal pressure (Mignon-Risse et al. 2021). Furthermore, on slightly larger scales, Peters et al. (2011) show that while magnetic fields reduce fragmentation, they do not prevent it and gravitational instability eventually sets in. Nevertheless, without information about the magnetic field strengths in the CORE sources, it is not trivial to estimate what extent that may affect our conclusions.

### 5.6.4. Disk thickness

One of the main assumptions for the Toomre analysis is the thin-disk approximation, such that $Q_{crit}$ for an isothermal disk of finite thickness is lower by ~0.7 (Goldreich & Lynden-Bell 1965). This is because the thicker a disk, the more diluted its self-gravity becomes. To get an estimate for the importance of disk thickness, Gammie (2001) derived an additional stability criterion taking into account the disk scale height, $H \approx c_s/\Omega$, such that a disk is unstable against gravitational instabilities if

$$M_{\text{disk}} \gtrsim \frac{H}{r} M_{*}. \qquad (13)$$

Accounting for the self-gravity of the disk (i.e. replacing $M_*$ with $M_* + M_{\text{disk}}$) and with the surface density of the disk $\Sigma = M_{\text{disk}}/(\pi r^2)$, the Toomre criterion can be rewritten as

$$Q = \left(\frac{H}{r}\right) \left(\frac{M_* + M_{\text{disk}}}{M_{\text{disk}}}\right). \qquad (14)$$

Using the disk scale height, $H \approx c_s/\Omega$, and Eqs. (10) and (12), we create maps of $H/r$ ratio and report the median $H/r$ values in Table 7. The observed median disk scale heights are generally too high (0.15–0.31) for a thin disk-like structure, possibly because at our resolution the envelope and disk components are likely blended, and the fact that temperatures derived from CH$_3$CN emission probe a layer higher than the disk midplane. Figure 11 shows the median Q as a function of $(H/r)(M_*+M_{\text{disk}})/M_{\text{disk}}$. We find a tight relationship between the two stability criteria, validating the analytically motivated scalings and the naïve expectations on the disk physics. However, we observe a normalisation value of ~2 (i.e. $Q \approx 2(H/r)(M_*+M_{\text{disk}})/M_{\text{disk}}$) instead of the expected value of unity, which is a direct consequence of the high disk scale heights.

## 6. Summary and conclusions

In this work, we present the gas kinematics for a sample of 20 high-mass star-forming regions observed with the NOEMA interferometer at a resolution of ~0.4" as part of the CORE survey. While the sample was chosen uniformly to target the early phase of high-mass star formation, there exists a diversity in age and chemical richness. We focussed our analysis on 15 sources for which the emission of the dense gas tracer CH$_3$CN was detected with the goal of searching for and characterising the properties of disk-like structures. Fig. 12 shows the main results that can be obtained for a given region with respect to the kinematics and derived physical properties. The following is a summary of our findings:

- Using the K-ladder rotational transitions of the dense gas tracer CH$_3$CN (12$_K$–11$_K$), we found velocity gradients across 13 cores perpendicular to the directions of bipolar molecular outflows, making them excellent disk candidates.
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Median *Q* plotted against gas mass (**left**) and stellar mass (**right**) for 13 candidate disks within the CORE survey, coloured according to the luminosity of the regions within which they reside. The Toomre-stable disks are marked by triangles.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Median *Q* as a function of ($\frac{H}{r}$)($\frac{M_{*}+M_{\text{disk}}}{M_{\text{disk}}}$). The Toomre-stable disks are marked by triangles.</figcaption></figure>

The directions of molecular outflows were determined from a combination of tracers ($^{12}$CO, $^{13}$CO, C$^{18}$O, and SO), both from single-dish IRAM 30-m observations as well as combined NOEMA and single-dish observations, included in the CORE survey.

- The rotating structures vary in size from ~1800 to ~8500 au, measured from the extent of CH$_3$CN (12$_3$–11$_3$) emission across the direction of the strongest velocity gradient.

- From the ratio of free-fall to rotational timescales as a function of gas mass, we established that the 13 sources in this sample are in fact rotationally supported disk candidates and not transient toroids.

- The velocity profiles of some disk candidates resemble differential rotation with high-velocity components closest to

the central object while others mimic rigid-body-like rotation. The latter may be a consequence of a lack of spatial resolution. We fitted the 5–7σ edge of the PV diagram for each disk with a Keplerian profile. The estimated dynamical masses, which for a Keplerian disk is the mass of the central object, are in the range ~10–25 $M_{\odot}$.

- Modelling the level population of CH$_3$CN (12$_K$–11$_K$) $K=0–6$ lines under LTE conditions using XCLASS, we provide temperature maps of the regions. Temperatures are on average warm (70–210 K, depending on the target) with a diversity in distributions. Some are warm in the centre and cooler in the outskirts, while others are more uniform in temperature. In some cases the ejection points of the outflows coincide with regions of warmer temperature where linewidths are also broader. This hints at the possibility that these regions have been heated by the outflows and associated shocks and/or regions that are carved out by the molecular outflows.

- We calculated the gas mass of the disk-like structures assuming gas and dust temperatures are in equilibrium and optimally thin dust emission. Gas masses range from 1 $M_{\odot}$ to as high as 30 $M_{\odot}$, depending on the region.

- Radial profiles of the specific angular momentum for the best disk candidates span a range of 1–2 orders of magnitude, on average around ~10⁻³ km s⁻¹ pc, and follow a $j \propto r^{1.7}$ relation. This scaling is consistent with synthetic observations of poorly resolved massive disk simulations where the envelope and disk components become completely blended. The radial profiles taper off at ~1000 au, suggesting that rotationally supported disks must exist below this scale. For the closest sources in our sample that are resolved at better linear resolutions, this limit was estimated at 300–400 au and 500–700 au for Cep A and IRAS 21078, respectively.

- Assuming the candidate disks to be in Keplerian rotation about a star with the estimated protostellar mass and accounting for the self-gravity of the disk by including its gas mass, we created maps of the Toomre *Q* parameter to determine their stability against axisymmetric gravitational collapse. Most disk candidates (11 out of 13) are Toomre unstable at the scales probed by our observations because of their high masses. In particular, disks with masses greater
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Summary plot showing the kinematics and derived properties for one source in the CORE sample (G75.78). **Top left panel:** Intensity-weighted mean velocity (first moment) map of CH₃CN (12₅–11₅) in colour with 1.37 mm continuum contours. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. The dotted line indicates the position of the strongest velocity gradient tracing the disk, i.e. perpendicular to the rotation axis. The edges of the assumed disk extent are marked with an ×. **Top right panel:** PV plot of CH₃CN (12₅–11₅) along the cut in the direction of rotation as depicted by the dotted line in the top left panel. The contours correspond to the 6σ level increasing in steps of 24σ. Yellow lines show the Keplerian rotation curves for an enclosed mass of 19 M⊙. The cross in the bottom right corner corresponds to the spatial and spectral resolutions. **Bottom left panel:** Rotational temperature maps obtained by fitting CH₃CN (12K–11K) K = 0–6 and CH₃¹³CN (12K–11K) K = 0–3 lines with XCLASS in colour and 1.37 mm continuum contours. **Bottom right panel:** Toomre Q map assuming a protostar is located at the position of the continuum peak as depicted by a star and accounting for the self-gravity of the disk, with 1.37 mm continuum contours. Regions outside of 6σ continuum contours are masked out. The synthesised beam and a scale bar are shown in the bottom of the panels.</figcaption></figure>

than ~10–20% of the mass of their host (proto)stars are subject to fragmentation through gravitational instabilities. We find that more luminous YSOs tend to have disks that are more massive and as a result more prone to fragmentation. These findings remain valid even when accounting for the thickness of the disks.

In conclusion, we find that disk-mediated accretion may be a common mechanism for high-mass star formation, similar to the formation scenario of their lower-mass counterparts. We further discover that most high-mass disk candidates in the CORE sample (11 out of 13) are prone to fragmentation due to gravitational instabilities early in their formation due to their high disk to stellar mass ratio. This has profound implications for the formation of the most massive stars, as it modifies the accretion evolution of high-mass protostars. Considering the fact that almost all high-mass stars are found in binary or multiple systems,

disk fragmentation seems to be an important mechanism in the formation of these systems.

**Acknowledgements.** The authors would like to thank the referee, Adam Ginsburg, whose comments helped improve the clarity and impact of the paper. A.A. and H.B. acknowledge support from the European Research Council under the European Community's Horizon 2020 framework programme (2014–2020) via the ERC Consolidator Grant 'From Cloud to Star Formation (CSF)' (project number 648505). RK 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. RGM acknowledges support from UNAM-PAPIIT project IN108822. AP acknowledges financial support from the UNAM-PAPIIT IN111421 grant, the Sistema Nacional de Investigadores of CONACYT. RGM and AP also acknowledge support from the CONACYT project number 86372 of the 'Ciencia de Frontera 2019' programme, entitled 'Citlal-cóatl: A multiscale study at the new frontier of the formation and early evolution of stars and planetary systems', México.
## References

Ahmadi, A., Beuther, H., Mottram, J. C., et al. 2018, A&A, 618, A46

Ahmadi, A., Kuiper, R., & Beuther, H. 2019, A&A, 632, A50

Ando, K., Nagayama, T., Omodaka, T., et al. 2011, PASJ, 63, 45

André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27

Apai, D., Bik, A., Kaper, L., Henning, T., & Zinnecker, H. 2007, ApJ, 655, 484

Bate, M. R. 2000, MNRAS, 314, 33

Bate, M. R., & Bonnell, I. A. 1997, MNRAS, 285, 33

Belloche, A. 2013, in EAS Publications Series, 62, 25

Beltrán, M. T., & de Wit, W. J. 2016, A&ARv, 24, 6

Beltrán, M. T., Cesaroni, R., Neri, R., & Codella, C. 2011, A&A, 525, A151

Beuther, H., Schilke, P., Gueth, F., et al. 2002, A&A, 387, 931

Beuther, H., Linz, H., & Henning, T. 2012, A&A, 543, A88

Beuther, H., Linz, H., & Henning, T. 2013, A&A, 558, A81

Beuther, H., Linz, H., Henning, T., Feng, S., & Teague, R. 2017, A&A, 605, A61

Beuther, H., Mottram, J. C., Ahmadi, A., et al. 2018, A&A, 617, A100

Boley, P. A., Linz, H., van Boekel, R., et al. 2013, A&A, 558, A24

Bosco, F., Beuther, H., Ahmadi, A., et al. 2019, A&A, 629, A10

Carter, M., Lazareff, B., Maier, D., et al. 2012, A&A, 538, A89

Cesaroni, R., Galli, D., Lodato, G., Walmsley, M., & Zhang, Q. 2006, Nature, 444, 703

Cesaroni, R., Sánchez-Monge, Á., Beltrán, M. T., et al. 2017, A&A, 602, A59

Cesaroni, R., Beuther, H., Ahmadi, A., et al. 2019, A&A, 627, A68

Chabrier, G. 2003, ASP, 115, 763

Chen, H.-R. V., Keto, E., Zhang, Q., et al. 2016, ApJ, 823, 125

Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925

Choi, Y. K., Hachisuka, K., Reid, M. J., et al. 2014, ApJ, 790, 99

Clark, B. G. 1980, A&A, 89, 377

Crowther, P. A. 2005, in IAU Symposium, 227, 389

Csenager, T., Bontemps, S., Wyrowski, F., et al. 2018, A&A, 617, A89

Davis, C. J., Moriarty-Schieven, G., Eislöffel, J., Hoare, M. G., & Ray, T. P. 1998, AJ, 115, 1118

De Buizer, J. M., Liu, M., Tan, J. C., et al. 2017, ApJ, 843, 33

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)

Elmegreen, B. G. 1991, ApJ, 378, 139

Elmegreen, B. G. 1994, ApJ, 433, 39

Fallscheer, C., Beuther, H., Zhang, Q., Keto, E., & Sridharan, T. K. 2009, A&A, 504, 127

Feng, S., Beuther, H., Semenov, D., et al. 2016, A&A, 593, A46

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, & T. Henning (Tucson: University of Arizona Press), 451

Gammie, C. F. 2001, ApJ, 553, 174

Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92

Gaume, R. A., Goss, W. M., Dickel, H. R., Wilson, T. L., & Johnston, K. J. 1995, ApJ, 438, 776

Gieser, C., Semenov, D., Beuther, H., et al. 2019, A&A, 631, A142

Gieser, C., Beuther, H., Semenov, D., et al. 2021, A&A, 648, A66

Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit Astrophysics Source Code Library, [record ascl:1109.0001]

Ginsburg, A., Koch, E., Robitaille, T., et al. 2019, https://doi.org/10.5281/zenodo.2573901

Ginsburg, A., Sokolov, V., de Val-Borro, M., et al. 2022, AJ, 163, 291

Goddi, C., Zhang, Q., & Moscadelli, L. 2015, A&A, 573, A108

Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125

Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

Guzmán, A. E., Sanhueza, P., Zapata, L., Garay, G., & Rodríguez, L. F. 2020, ApJ, 904, 77

Hacar, A., Clark, S., Heitsch, F., et al. 2022, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura

Hachisuka, K., Brunthaler, A., Menten, K. M., et al. 2006, ApJ, 645, 337

Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25

Henning, T., Lapinov, A., Schreyer, K., Stecklum, B., & Zinchenko, I. 2000, A&A, 364, 613

Ilee, J. D., Cyganowski, C. J., Nazari, P., et al. 2016, MNRAS, 462, 4386

Ilee, J. D., Cyganowski, C. J., Brogan, C. L., et al. 2018, ApJ, 869, L24

Jiménez-Serra, I., Martín-Pintado, J., Caselli, P., et al. 2009, ApJ, 703, L157

Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19

Johnston, K. G., Hoare, M. G., Beuther, H., et al. 2020, A&A, 634, A11

Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., I., & Lee, C. W. 2008, A&A, 487, 993

Kennicutt, Robert C., J. 1989, ApJ, 344, 685

Keto, E. R. 1991, ApJ, 371, 163

Kim, W.-T., & Ostriker, E. C. 2001, ApJ, 559, 70

Klaassen, P. D., Wilson, C. D., Keto, E. R., et al. 2011, A&A, 530, A53

Klassen, M., Pudritz, R. E., Kuiper, R., Peters, T., & Banerjee, R. 2016, ApJ, 823, 28

Kölligan, A., & Kuiper, R. 2018, A&A, 620, A182

Kratter, K. M., & Lodato, G. 2016, ARA&A, 54, 271

Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563

Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585

Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005, ApJ, 618, L33

Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754

Kuiper, R., & Hosokawa, T. 2018, A&A, 616, A101

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20

Kuiper, R., Yorke, H. W., & Turner, N. J. 2015, ApJ, 800, 86

Kuiper, R., Turner, N. J., & Yorke, H. W. 2016, ApJ, 832, 40

Leurini, S., Codella, C., Zapata, L., et al. 2011, A&A, 530, A12

Lumsden, S. L., Hoare, M. G., Urquhart, J. S., et al. 2013, ApJS, 208, 11

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Maud, L. T., Moore, T. J. T., Lumsden, S. L., et al. 2015, MNRAS, 453, 645

Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2018, A&A, 620, A31

Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2019, A&A, 627, A6

Mège, P., Russeil, D., Zavagno, A., et al. 2021, A&A, 646, A74

Mestel, L. 1963, MNRAS, 126, 553

Meyer, D. M. A., Kuiper, R., Kley, W., Johnston, K. G., & Vorobyov, E. 2018, MNRAS, 473, 3615

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2021, A&A, 652, A69

Molinari, S., Brand, J., Cesaroni, R., & Palla, F. 1996, A&A, 308, 573

Molinari, S., Testi, L., Brand, J., Cesaroni, R., & Palla, F. 1998, ApJ, 505, L39

Möller, T., Endres, C., & Schilke, P. 2017, A&A, 598, A7

Moscadelli, L., & Goddi, C. 2014, A&A, 566, A150

Moscadelli, L., Reid, M. J., Menten, K. M., et al. 2009, ApJ, 693, 406

Moscadelli, L., Sanna, A., Cesaroni, R., et al. 2019, A&A, 622, A206

Moscadelli, L., Beuther, H., Ahmadi, A., et al. 2021, A&A, 647, A114

Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243

Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41

Mottram, J. C., Hoare, M. G., Urquhart, J. S., et al. 2011, A&A, 525, A149

Mottram, J. C., Beuther, H., Ahmadi, A., et al. 2020, A&A, 636, A118

Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2022, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, and M. Tamura

Ohashi, N., Hayashi, M., Ho, P. T. P., & Momose, M. 1997a, ApJ, 475, 211

Ohashi, N., Hayashi, M., Ho, P. T. P., et al. 1997b, ApJ, 488, 317

Oliva, G. A., & Kuiper, R. 2020, A&A, 644, A41

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Peters, T., Banerjee, R., Klessen, R. S., et al. 2010a, ApJ, 711, 1017

Peters, T., Klessen, R. S., Mac Low, M.-M., & Banerjee, R. 2010b, ApJ, 725, 134

Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72

Pineda, J. E., Zhao, B., Schmiedeke, A., et al. 2019, ApJ, 882, 103

Pomohaci, R., Oudmaijer, R. D., & Goodwin, S. P. 2019, MNRAS, 484, 226

Pudritz, R. E., & Ray, T. P. 2019, Front. Astron. Space Sci., 6, 54

Qiu, K., Zhang, Q., & Menten, K. M. 2011, ApJ, 728, 6

Ramírez-Tannus, M. C., Backs, F., de Koter, A., et al. 2021, A&A, 645, A10

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131

Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79

Sakai, N., Reid, M. J., Menten, K. M., Brunthaler, A., & Dame, T. M. 2019, ApJ, 876, 30

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444

Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15

Sánchez-Monge, Á., Beltrán, M. T., Cesaroni, R., et al. 2014, A&A, 569, A11

Sandell, G., Goss, W. M., Wright, M., & Corder, S. 2009, ApJ, 699, L31

Sanna, A., Moscadelli, L., Surcis, G., et al. 2017, A&A, 603, A94

Sanna, A., Kölligan, A., Moscadelli, L., et al. 2019, A&A, 623, A77

Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320

Seifried, D., Sánchez-Monge, Á., Walch, S., & Banerjee, R. 2016, MNRAS, 459, 1892

Steer, D. G., Dewdney, P. E., & Ito, M. R. 1984, A&A, 137, 159

Suri, S., Beuther, H., Gieser, C., et al. 2021, A&A, 655, A84

Tanaka, K. E. I., Zhang, Y., Hirota, T., et al. 2020, ApJ, 900, L2

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Tobin, J. J., Hartmann, L., Bergin, E., et al. 2012, ApJ, 748, 16

Toomre, A. 1964, ApJ, 139, 1217

Urquhart, J. S., Morgan, L. K., Figura, C. C., et al. 2011, MNRAS, 418, 1689
Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Xu, Y., Reid, M. J., Zheng, X. W., & Menten, K. M. 2006, Science, 311, 54

Xu, Y., Reid, M. J., Menten, K. M., et al. 2009, ApJ, 693, 413

Xu, Y., Li, J. J., Reid, M. J., et al. 2013, ApJ, 769, 15

Xu, Y., Reid, M., Dame, T., et al. 2016, Sci. Adv., 2, e1600878

Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846

Zapata, L. A., Palau, A., Galván-Madrid, R., et al. 2015, MNRAS, 447, 1826

Zapata, L. A., Garay, G., Palau, A., et al. 2019, ApJ, 872, 176

Zhang, Q., Hunter, T. R., Brand, J., et al. 2001, ApJ, 552, L167

Zhang, Q., Hunter, T. R., Brand, J., et al. 2005, ApJ, 625, 864

Zhang, S., Hartmann, L., Zamora-Avilés, M., & Kuznetsova, A. 2018, MNRAS, 480, 5495

Zinchenko, I., Liu, S.-Y., Su, Y.-N., & Wang, Y. 2018, Res. Astron. Astrophys., 18, 093

¹ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
e-mail: aahmadi@strw.leidenuniv.nl

² Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

³ Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1, 85749 Garching bei München, Germany

⁴ Department of Astrophysics, University of Vienna, Türkenschanzstraße 17 (Sternwarte) 1180 Wien, Austria

⁵ Faculty of Physics, University of Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany

⁶ Institut de Ciències de l'Espai (ICE, CSIC), Can Magrans s/n, 08193, Bellaterra, Barcelona, Spain

⁷ Institut d'Estudis Espacials de Catalunya (IEEC), Barcelona, Spain

⁸ Observatorio Astronómico Nacional (OAN, IGN), Calle Alfonso XII 3, 28014, Madrid, Spain

⁹ I. Physikalisches Institut der Universität zu Köln, Zülpicher Str. 77, 50937, Köln, Germany

¹⁰ Dept. of Physics and Astronomy, McMaster University, 1280 Main Street West, Hamilton, ON L8S 4K1, Canada

¹¹ Department of Chemistry, Ludwig Maximilian University, Bütendandstr. 5-13, 81377 Munich, Germany

¹² IRAM, 300 rue de la Piscine, Domaine Universitaire, 38406 Saint-Martin-d'Hères, France

¹³ INAE – Osservatorio Astrofísico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

¹⁴ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France

¹⁵ Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Antigua Carretera a Pátzcuaro 8701, Ex-Hda. San José de la Huerta, 58089, Morelia, Michoacán, Mexico

¹⁶ School of Physics and Astronomy, The University of Leeds, Sir William Henry Bragg Building, Woodhouse Lane, Leeds, LS2 9JT, UK

¹⁷ Center for Astrophysics, Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02420, USA

¹⁸ UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK

¹⁹ INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047, Selargius (CA), Italy

²⁰ Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK

²¹ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748, Garching bei München, Germany

²² Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany

²³ School of Physics and Astronomy, Cardiff University, Queen's Buildings, The Parade, Cardiff, CF24 3AA, UK

²⁴ Centre for Astrophysics and Planetary Science, University of Kent, Canterbury, CT2 7NH, UK

²⁵ Universidad Autónoma de Chile, Nucleo de Astroquímica y Astrofísica, 425 Avda Pedro de Valdivia, Providencia, Santiago de Chile, Chile
## Appendix A: Moment maps

In this appendix, we present zeroth and second moment maps of CH₃CN (12₃ – 11₃) in figures A.1 and A.2, respectively.

<figure><img src="image_15.png" /><figcaption>Fig. A.1: Integrated intensity (zeroth moment) maps of CH₃CN (12₃ – 11₃) showing the dense gas distribution for 15 of the 20 sources in the CORE survey. The contours and features are as described in Fig. 2.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. A.2: Intensity-weighted velocity dispersion (second moment) maps of CH$_3$CN (12$_3$ − 11$_3$) showing the dense gas kinematics for 15 of the 20 sources in the CORE survey. The contours and features are as described in Fig. 2.</figcaption></figure>
## Appendix B: Decomposed velocity maps

We have fitted Gaussian profiles to the spectra of all sources using the *Spectral-Cube* tool within the *PySpecKit* package in Python (Ginsburg & Mirocha 2011; Ginsburg et al. 2019) and present the decomposed peak intensity, peak velocity, and full-width-at-half-maximum (FWHM) linewidth maps in this section (Figs. B.1 – B.3). Fitting one Gaussian profile assumes that there is only one velocity component at each position.

<figure><img src="image_17.png" /></figure>

Fig. B.1: Peak intensity (amplitude) maps of CH₃CN (12₃⁻¹¹₃) showing the dense gas kinematics for 15 of the 20 sources in the CORE survey obtained by fitting Gaussian profiles to their spectra. The contours and features are as described in Fig. 2.
<figure><img src="image_18.png" /></figure>

Fig. B.2: Peak velocity maps of CH$_3$CN (12$_3$ – 11$_3$) showing the dense gas kinematics for 15 of the 20 sources in the CORE survey obtained by fitting Gaussian profiles to their spectra. The contours and features are as described in Fig. 2.
<figure><img src="image_19.png" /></figure>

Fig. B.3: Linewidth (FWHM) maps of CH$_3$CN (12$_3$ − 11$_3$) showing the dense gas kinematics for 15 of the 20 sources in the CORE survey obtained by fitting Gaussian profiles to their spectra. The contours and features are as described in Fig. 2.
Appendix C: Molecular outflows

In this section, the distribution of outflow and shock tracing transitions are presented in Figs C.1 – C.6. The absolute difference between the assumed disk and outflow position angles is shown in Fig. C.7.

<figure><img src="image_20.png" /></figure>

Fig. C.1: Intensity maps of CO (2–1) emission from IRAM 30-m telescope integrated over the blueshifted and redshifted wings of emission, showing the outflow structure. The position of the strongest source in the field is depicted by a star. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. The map size of the IRAM 30-m observations is 1.5' by 1.5' with a half-power beam width of ~11" at this frequency.
<figure><img src="image_21.png" /></figure>

Fig. C.2: Intensity maps of ¹³CO (2 – 1) emission from IRAM 30-m telescope integrated over the blueshifted and redshifted wings of emission, showing the outflow structure. The position of the strongest source in the field is depicted by a star. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. The map size of the IRAM 30-m observations is 1.5' by 1.5' with a half-power beam width of ~11'' at this frequency.
<figure><img src="image_22.png" /></figure>

Fig. C.3: Intensity maps of ¹³CO (2–1) emission from merged NOEMA and IRAM 30-m data integrated over the blueshifted and redshifted wings of emission, showing the outflow structure. The position of the strongest source in the field is depicted by a star. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. We did not merge the data for the pilot sources as we did not observe this target with the most compact NOEMA configuration (D-array), hence the missing panels.
<figure><img src="image_23.png" /></figure>

Fig. C.4: NOEMA 1.37 mm continuum images (greyscale) and intensity maps of ¹³CO (2 – 1) integrated over the blueshifted and redshifted wings of emission (blue and red contours), tracing either outflows or disk winds. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. Note that most of the emission is filtered out by the interferometer.
<figure><img src="image_24.png" /></figure>

Fig. C.5: Intensity maps of C$^{18}$O (2 – 1) emission from IRAM 30-m telescope integrated over the blueshifted and redshifted wings of emission, showing the outflow structure. The position of the strongest source in the field is depicted by a star. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. The map size of the IRAM 30-m observations is 1.5' by 1.5' with a half-power beam width of ~11'' at this frequency. C$^{18}$O is not a good outflow tracer for all sources, hence the lack of contours in some targets.
<figure><img src="image_25.png" /></figure>

Fig. C.6: Intensity maps of SO(6<sub>5</sub> − 5<sub>4</sub>) emission from IRAM 30-m telescope integrated over the blueshifted and redshifted wings of emission, showing the outflow structure. The position of the strongest source in the field is depicted by a star. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. A scale-bar is shown in the bottom right corner of each panel. The map size of the IRAM 30-m observations is 1.5' by 1.5' with a half-power beam width of ~11'' at this frequency.
<figure><img src="image_26.png" /><figcaption>Fig. C.7: Distribution of the absolute difference between the disk position angles and the assumed outflow position angles for our sample of 13 disk candidates in the CORE survey.</figcaption></figure>

# Appendix D: NGC7538 IRS1

## Appendix D.1: Rotational temperature

In NGC7538 IRS1, the molecular lines are seen in absorption against the strong continuum source with many velocity components due to the complex dynamics in this source. To properly model the temperature structure across this disk candidate, we used XCLASS to model the molecular line emission existing outside the continuum to the south-west where the spectra are not affected by absorption features. We then scaled the modelled rotational temperature from the edge of the continuum structure (~135 K) towards the continuum peak following the temperature power-law distribution $T \propto r^{-0.4}$, derived by Gieser et al. (2021) from the CORE survey. The resulting rotational temperature map is shown in Fig. 3.

## Appendix D.2: Free-free contribution

NGC7538 IRS1 is the only source in our sample of disk candidates with significant free-free contribution in the millimetre (Beuther et al. 2017). For this source, we smoothed the continuum emission map of Beuther et al. (2017) observed with the Karl G. Jansky Very Large Array (VLA) at 24.7 GHz (1.2 cm) to the resolution of our observations. We then scaled the emission of this smoothed map by $y^{-0.1}$, the expected frequency-dependence of free-free emission, and subtracted it from our millimetre continuum observations. The continuum flux densities before and after the subtraction are provided in the body of Table 4 and as a note under the table, respectively. The contribution is ~13%.

# Appendix E: XCLASS maps

XCLASS output maps such as the column density, velocity offset, linewidth, and source size maps for the best disk candidates in the CORE survey are shown in Figs. E.1–E.4. The XCLASS output maps for the remaining three sources that had CH$_3$CN emission but that were not good disk candidates are presented in Fig. E.5.
<figure><img src="image_27.png" /></figure>

Fig. E.1: CH₃CN column density maps obtained by fitting CH₃CN (12<sub>K</sub> - 11<sub>K</sub>) K = 0–6 and CH₃<sup>13</sup>CN (12<sub>K</sub> - 11<sub>K</sub>) K = 0–3 lines with XCLASS for the best disk candidates in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. For NGC7538 IRS1, only the region outside the continuum to the south-west is modelled by XCLASS as the molecular lines are seen in absorption against the strong continuum source with many velocity components due to the complex dynamics in this source. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.
<figure><img src="image_28.png" /></figure>

Fig. E.2: Maps of offset velocity with respect to the systemic velocity obtained by fitting CH₃CN(12K – 11K) K = 0–6 and CH₃¹³CN(12K – 11K) K = 0–3 lines with XCLASS for the best disk candidates in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. For NGC7538 IRS1, only the region outside the continuum to the south-west is modelled by XCLASS as the molecular lines are seen in absorption against the strong continuum source with many velocity components due to the complex dynamics in this source. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.
<figure><img src="image_29.png" /></figure>

Fig. E.3: Maps of linewidth obtained by fitting CH₃CN (12K – 11K) K = 0 – 6 and CH₃¹³CN (12K – 11K) K = 0 – 3 lines with XCLASS for the best disk candidates in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. For NGC7538 IRS1, only the region outside the continuum to the south-west is modelled by XCLASS as the molecular lines are seen in absorption against the strong continuum source with many velocity components due to the complex dynamics in this source. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.
<figure><img src="image_30.png" /></figure>

Fig. E.4: Maps of source size parameter obtained by fitting CH₃CN (12<sub>K</sub> - 11<sub>K</sub>) K = 0-6 and CH₃<sup>13</sup>CN (12<sub>K</sub> - 11<sub>K</sub>) K = 0-3 lines with XCLASS for the best disk candidates in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. For NGC7538 IRS1, only the region outside the continuum to the south-west is modelled by XCLASS as the molecular lines are seen in absorption against the strong continuum source with many velocity components due to the complex dynamics in this source. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.
<figure><img src="image_31.png" /></figure>

Fig. E.5: Maps of rotational temperature (top row), CH₃CN column density (second row), offset velocity (third row), linewidth (fourth row), and source size parameter (bottom row) obtained by fitting CH₃CN (12K – 11K) K = 0 – 6 and CH₃¹³CN (12K – 11K) K = 0 – 3 lines with XCLASS for the remaining three sources (columns) in the CORE survey. The contours correspond to the 1.37 mm continuum as described in Fig. 1. The blue and red arrows correspond to the estimated directions of bipolar blueshifted and redshifted molecular outflows, respectively. The synthesised beam is shown in the bottom left corner and a scale bar in the bottom right corner of each panel.
## Appendix F: Specific angular momentum

In order to better understand the observational findings, we analyse the expected specific angular momentum profiles for the numerical simulations presented in Oliva & Kuiper (2020) to study the effects of inclination and angular resolution on the specific angular momentum, *j*. The 3D radiation-hydrodynamic simulation under investigation follows the collapse of a massive core with 200 $M_\odot$ in gas and dust, which results in the formation and fragmentation of a massive disk. At this snapshot, the star has gained 10 $M_\odot$, being fed from a fragmented disk with a mass of 8 $M_\odot$. It is important to note that the simulation presented here does not include magnetic fields, which are an important ingredient for disk formation. It also does not include an inhomogeneous, large-scale accretion flow from which the disk is fed, and an initial axially symmetric rotation profile of $\Omega(r) \propto r^{-0.75}$ was adopted.

We presented synthetic observations of this system in Ahmadi et al. (2019) and here calculate the observed $j \sin(i)$ according to Eq. (7) using the first moment maps of CH$_3$CN (12$_3$ – 11$_3$) at an inclination angle of 10°, 30°, 60°, and 80° from the plane of the sky along a cut in the E-W direction for the redshifted and blueshifted sides of the disk (see Fig. 5 of Ahmadi et al. 2019). Figure F.1 shows the local specific angular momentum radial profiles for the simulations at the different inclinations as well as for a set of synthetic observations with successively lower resolution (as noted in the panels). Note that we have purposefully not corrected for the known inclination in order to compare the findings with our NOEMA observations. Similar to the observational findings for low-mass star formation, we find three distinct regimes in the *j* radial profile. Least-squared power-law fit parameters to these three distinct regimes in the specific angular momentum plots for the numerical simulations are listed in Table F.1. The inner disk region ($r < 200$ au) actively accreting material onto the central protostar with $j \propto r^{0.6}$ close to Keplerian rotation ($j \propto r^{0.5}$), an intermediate region where *j* is quite flat but not exactly constant, and the region beyond 700 au is rather consistent with the $j \propto r^{1.6}$ relation that had been found for the rotation of low-mass dense cores and Class 0 envelopes (e.g. Gaudel et al. 2020). The intermediate regime has a significant amount of substructure as four fragments reside in this range of radii at this snapshot, accreting part of the disk material. The discontinuity at ~700 au roughly corresponds to the position where the rotating core or envelope begins to dynamically collapse. We find that the simulation with more of a face-on view has almost an order of magnitude lower specific angular momentum with flatter radial profiles and for which the distinction between the different regimes is unclear.

As the angular resolution worsens, the substructure in the *j* radial profile becomes smoother with the power-law relations for the different regimes becoming more continuous and not as clear to disentangle. The specific angular momentum of the blueshifted side seems to have a more pronounced kink than the redshifted side, best visible in the synthetic observations at 800 pc. This is due to the existence of a fragment at this position, accreting some of the disk material and complicating the angular momentum profile at this position. At the scale of synthetic NOEMA observations at 2 kpc, the envelope and disk components become completely blended, with only one power-law relation seen across all scales, with $j \propto r^{1.6}$. Interestingly, there still exists a small shallowing of the slope (a kink) in the *j* profiles.

<table><caption>Table F.1: Fit parameters for the different regions of the specific angular momentum radial profiles of the simulations shown in the top left panel of Fig. F.1.</caption><thead><tr><th rowspan="2">Inclination</th><th colspan="4">Inner (r &lt; 100 au)</th></tr><tr><th>a<sub>blue</sub></th><th>a<sub>red</sub></th><th>b<sub>blue</sub></th><th>b<sub>red</sub></th></tr></thead><tbody><tr><td>10°</td><td>0.56 ± 0.06</td><td>0.56 ± 0.06</td><td>-3.16 ± 0.06</td><td>-3.16 ± 0.06</td></tr><tr><td>30°</td><td>0.54 ± 0.07</td><td>0.54 ± 0.07</td><td>-2.71 ± 0.07</td><td>-2.71 ± 0.07</td></tr><tr><td>60°</td><td>0.56 ± 0.09</td><td>0.56 ± 0.09</td><td>-2.6 ± 0.1</td><td>-2.6 ± 0.1</td></tr><tr><td>80°</td><td>0.75 ± 0.08</td><td>0.75 ± 0.08</td><td>-2.51 ± 0.08</td><td>-2.51 ± 0.08</td></tr></tbody><thead><tr><th rowspan="2">Inclination</th><th colspan="4">Middle (100 au &lt; r &lt; 400 au)</th></tr><tr><th>a<sub>blue</sub></th><th>a<sub>red</sub></th><th>b<sub>blue</sub></th><th>b<sub>red</sub></th></tr></thead><tbody><tr><td>10°</td><td>0.14 ± 0.05</td><td>-0.01 ± 0.08</td><td>-3.36 ± 0.03</td><td>-3.77 ± 0.05</td></tr><tr><td>30°</td><td>0.15 ± 0.04</td><td>0.26 ± 0.04</td><td>-2.96 ± 0.02</td><td>-3.04 ± 0.02</td></tr><tr><td>60°</td><td>0.35 ± 0.03</td><td>0.54 ± 0.05</td><td>-2.69 ± 0.02</td><td>-2.7 ± 0.03</td></tr><tr><td>80°</td><td>0.52 ± 0.04</td><td>0.48 ± 0.03</td><td>-2.58 ± 0.02</td><td>-2.83 ± 0.02</td></tr></tbody><thead><tr><th rowspan="2">Inclination</th><th colspan="4">Outer (r &gt; 400 au)</th></tr><tr><th>a<sub>blue</sub></th><th>a<sub>red</sub></th><th>b<sub>blue</sub></th><th>b<sub>red</sub></th></tr></thead><tbody><tr><td>10°</td><td>0.78 ± 0.03</td><td>0.6 ± 0.1</td><td>-3.12 ± 0.01</td><td>-3.53 ± 0.02</td></tr><tr><td>30°</td><td>0.74 ± 0.05</td><td>1.0 ± 0.1</td><td>-2.76 ± 0.01</td><td>-2.85 ± 0.02</td></tr><tr><td>60°</td><td>1.30 ± 0.05</td><td>1.75 ± 0.09</td><td>-2.45 ± 0.01</td><td>-2.45 ± 0.01</td></tr><tr><td>80°</td><td>1.45 ± 0.03</td><td>1.75 ± 0.05</td><td>-2.34 ± 0.01</td><td>-2.34 ± 0.01</td></tr></tbody></table>

Notes. The heading parameters *a* and *b* satisfy the following relation: $j \sin(i) = 10^b (\frac{r}{1000 \text{ au}})^a \text{ km s}^{-1} \text{ pc}$, for the blueshifted and redshifted sides.
<figure><img src="image_32.png" /></figure>

Fig. F.1: Specific angular momentum radial profiles (see Eq. 7) for a high-resolution numerical simulation of a high-mass core having formed a massive disk that is feeding the central protostar, inclined to 10°, 30°, 60°, and 80° (top left), their synthetic ALMA observations at 2000 pc (top right), and synthetic NOEMA observations at 800 pc (bottom left) and 2000 pc (bottom right), (see Ahmadi et al. 2019). The blue and red colours correspond to the blueshifted and redshifted sides, respectively, and the hues correspond to the different inclinations labelled in the bottom-right panel. Note that we have not corrected for the known inclination in order to compare the simulated observations to our NOEMA observations with unknown inclinations. The distance between each neighbouring point is a half beam spacing. The grey dashed lines show various slopes to guide the eye.