© ESO 2021

1 Introduction

High-mass (≳8 M_⊙ and ≳10⁴ L_⊙) stars form within young stellar object (YSO) clusters deeply embedded inside thick dust cocoons, whose study requires the superior angular resolution and sensitivity achieved only recently by the new or upgraded (sub)millimeter interferometers, such as the Atacama Large Millimeter/submillimeter Array (ALMA), the NOrthern Extended Millimeter Array (NOEMA), and the SubMillimeter Array (SMA). The processes of the formation of these clusters via the collapse and fragmentation of the parental molecular clump into multiple smaller cores, and mass accretion on individual cores, are not yet clear. Observing similar scales of 0.1–1 pc with comparable angular resolution of ≳1′′, recent SMA and ALMA studies of infrared-quiet massive clumps have revealed different fragmentation properties: in some cases, limited fragmentation with a large fraction of cores with masses in the range 8–120 M_⊙ (Wang et al. 2014; Csengeri et al. 2017; Neupane et al. 2020), and, in other cases, a large population of low-mass (≤1 M_⊙) cores with a maximums core mass of 11 M_⊙ (Sanhueza et al. 2019). There is a clear need for extending this kind of studies to identify the main agent(s) of molecular clump fragmentation. In the very early evolutionary phases, one would expect negligible support against gravitational collapse from thermal pressure and protostellar feedback (such as outflows or HII regions), and theoretical models predict that the main competitors of gravity should be the internal clump turbulence and magnetic pressure (e.g., Federrath 2015; Klessen & Glover 2016; Hennebelle et al. 2019).

The fragmentation modality determines the mass reservoir for the formation of individual stars and affects the accretion process as well, as discussed by the two main competing theories: the “core-accretion” model (McKee & Tan 2003), which predicts the existence of quasi-equilibrium massive cores providing all the material to form high-mass stars, and the “competitive accretion” model (Bonnell et al. 2004), in which the parent clump fragments into many low-mass cores that competitively accrete the surrounding gas. Within a core cluster, tidal interactions among nearby (~1000 au) YSOs can have a strong impact on accretion disks, causing them to precess (Cesaroni et al. 2005) and to suffer perturbations (Winter et al. 2018) or truncation (Goddi et al. 2011), depending on the relative mass and separation of the interacting YSOs.

Surveys at infrared and radio wavelengths have discovered that filaments are ubiquitous in massive molecular clouds (Molinari et al. 2010a; Lu et al. 2018) and have a wide range in lengths (a few to 100 pc) and line masses (a few hundreds to thousands of M_⊙ pc⁻¹). Young stellar object clusters and high-mass YSOs are often found at the junctions (named “hubs”) of this filamentary structures in molecular clouds (see, for instance, Myers 2009), and longitudinal mass flows (with typical rates of 10⁻⁴ –10⁻³ M_⊙ yr⁻¹) along the filaments converging toward the hubs have been identified (Chen et al. 2019; Schwörer et al. 2019; Treviño-Morales et al. 2019). These recent findings strongly suggest that hub-filament systems can play a fundamental role in the formation of high-mass stars and have led to complementary or alternative theories to the core-accretion and competitive-accretion models, such as the “global hierarchical collapse” (Vázquez-Semadeni et al. 2019), “conveyor belt” (Longmore et al. 2014), and “inertial-inflow” (Padoan et al. 2020) models. The main difference among these theories regards the origin of the hub-filament structures inside molecular clouds and the main driver (turbulence, density gradient, cloud-cloud collision, etc.) of the mass flows, but they all agree that the material to form the YSO clusters can be gathered from very large scales of 1–10 pc, channeled through the molecular cloud filaments.

Aiming at a statistical study of clump fragmentation and disk properties in high-mass YSOs, we are carrying out the large program “CORE” (P.I.: Henrik Beuther; team web-page¹), observing a sample of 20high-mass star-forming regions with the IRAM interferometer NOEMA in the 1.37 mm continuum and line emissions at high-angular resolution (0.′′4). The analysis of the continuum emission from cold dust shows a large variety in the fragmentation properties of the sample, from regions where a single massive core dominates to regions where as many as 20 cores, with comparable masses, are detected(Beuther et al. 2018). In several CORE targets, employing dense-gas molecular tracers such as the CH₃ CN and CH₃ OH rotational transitions, LSR velocity (V_LSR) gradients are often detected toward the most massive cores of the cluster (Ahmadi et al. 2018; Bosco et al. 2019; Cesaroni et al. 2019; Gieser et al. 2019). Extending over typical linear scales of ≲1000 au, these V_LSR gradients are interpreted in terms of the edge-on rotation of either a disk around a high-mass YSO or an unresolved binary (or multiple) system, or the combination of both motions.

IRAS 21078+5211, at a distance of 1.63 ± 0.05 kpc (Xu et al. 2013), is one of the closest and, with a bolometric luminosity of 1.3 × 10⁴ L_⊙, least luminous CORE sources (Beuther et al. 2018, see Table 1). IRAS 21078+5211 (alias G092.69+3.08) also belongs to the target list of the Protostellar Outflows at the EarliesT Stages (POETS) survey (Moscadelli et al. 2016; Sanna et al. 2018), which has recently been carried out with the goal of imaging the inner (10–100 au) outflow scales in a sample of 38 high-mass YSOs. Employing Very Long Baseline Array (VLBA) observations of the 22 GHz water masers and Jansky Very Large Array (JVLA) multi-frequency continuum observations, we investigated both the molecular and ionized component of the outflows. In IRAS 21078+5211, a jet is detected in the radio continuum emission, consisting of two slightly resolved components separated by ≈ 0.′′5 along the SW-NE direction (Moscadelli et al. 2016, see Fig. 8). While the NE component has a positive spectral index and traces the ionized core of the jet close to the YSO, the SW lobe has a negative spectral index owing to synchrotron emission. The 22 GHz water masers are located to the NE from the YSO at a separation between 100 and 300 au. The masers are clearly tracing the jet kinematics because their positions and proper motions are collimated along the radio jet.

In this paper we combine the CORE and POETS data to study the process of star formation (SF) in IRAS 21078+5211 over linear scales from 10⁴ to 10 au. Recent Large Binocular Telescope (LBT) observations are also employed to reveal the structure of the molecular outflows at the largest scales. The CORE and LBT observational setups, data reduction and analysis are described in Sect. 2, including a brief summary of the POETS observations of IRAS 21078+5211. In Sect. 3, we report on the main observational results going from the large scale, that is the distribution and kinematics of the molecular cores in the cluster, to the small scale, that is the detailed view of the accretion-ejection structures around the most massive YSO(s). In Sect. 4, these results are interpreted within the frame of a coherent picture connecting the multi-scale physical processes. Finally, our conclusions are drawn in Sect. 5.

2 Observations

2.1 CORE

The CORE program employs IRAM/NOEMA multi-configuration interferometric observations complemented with short-spacing IRAM 30 m single-dish data. A full description of the IRAM/NOEMA observational strategy and sample selection are given in Beuther et al. (2018) and the IRAM 30 m observations are detailed in Ahmadi et al. (2018) and Mottram et al. (2020). Hereafter, we report on the observations, data reduction, and analysis specifically for IRAS 21078+5211.

2.1.1 NOEMA

The IRAM/NOEMA interferometer observed IRAS 21078+5211 at 1.37 mm, in track sharing with another CORE target G100.3779−03.578, in three configurations (A, B, and D) between December 2014 and December 2016. The phase center for IRAS 21078+5211 was RA(J2000): 21^h 09^m 21.s64 and Dec(J2000): +52° 22′ 37.′′50, and the systemic V_LSR was assumed equal to V_sys = −6.1 km s⁻¹. In the course of the CORE observing campaign, the number of the NOEMA antennas increased from 6 to 8. The baseline lengths in the uv -plane range from 34 to 765 m, corresponding to spatial frequencies from 0.′′37 to 8.′′ 3. The gain calibrators were the quasars J2201+508 and 2037+511, and 3C454.3 and MWC349 were the bandpass and flux calibrators.

By employing the WideX correlator, we recorded a broad band covering the frequency range 217.167–220.834 GHz at 1.95 MHz spectral resolution (corresponding to ≈ 2.7 km s⁻¹), and eight narrow-band high-spectral resolution (≈ 0.43 km s⁻¹) units distributed over the broad band. The wide band is used to extract the line-free continuum, to get a chemical census of the region, and to study the distribution of more diffuse gas and the outflow kinematics (through low-excitation lines of the ¹³CO, H₂ CO and SO molecules). The eight narrow-band spectral units are centered at specific spectral locations to observe typical dense-gas tracers (such as high-excitation lines of CH₃CN, CH₃ OH, and CH₂ CO) for investigating the gas kinematics and physical conditions at high-angular resolution. A detailed description of the spectral line coverage can be found in Beuther et al. (2018, see Table 8) and Ahmadi et al. (2018). Table 1 shows the parameters, taken from the CDMS (Cologne Database for Molecular Spectroscopy², Müller et al. 2001, 2005) and JPL (Jet Propulsion Laboratory Catalog of Molecular Spectroscopy³, Pickett et al. 1998) databases, of all the molecular lines analyzed in this article.

The NOEMA data were calibrated with the CLIC and imaged with the MAPPING package in gildas⁴. The continuum data were extracted from line-free channels of the wide-band WideX spectra. Phase self-calibration was performed on the continuum data using the SELFCAL procedure. The gain table containing the self-calibration solutions was then applied to the narrow- and wide-band spectral line data using the UV_CAL task. A detailed description of the phase self-calibration of the CORE continuum and spectral line data can be found in Gieser et al. (2021).

The NOEMA continuum data were imaged with uniform weighting using the Clark algorithm (Clark 1980). The synthesized beam of the continuum image of IRAS 21078+5211 has full width at half maximum (FWHM) major and minor sizes of 0.′′48 and 0.′′ 33, and PA = 41°. For the extended ¹³CO and SO emissions, we used the low-resolution WideX spectral line data smoothed to a spectral resolution of 3.0 km s⁻¹. In order to avoid missing flux due to missing short-spacings, the data were combined with the IRAM 30 m observationsand imaged with a robust parameter of 3 using the SDI algorithm (Steer et al. 1984). For IRAS 21078+5211, the synthesized beam of the images produced by merging the NOEMA and IRAM 30 m spectral line data has FWHM major and minor sizes of 0.′′70 and 0.′′ 61, and PA = 66°. For the compact emissions of CH₃CN, HC₃N, and CH₂CO, we used the narrow-band NOEMA-only data smoothed to a spectral resolution of 0.5 km s⁻¹ and imaged with uniform weighting using the Clark algorithm. The NOEMA images of IRAS 21078+5211 have FWHM major and minor beam sizes of 0.′′46 and 0.′′ 31, with beam PA = 40°.

2.1.2 IRAM 30 m

IRAS 21078+5211 was observed with the IRAM 30 m telescope in the 1.4 mm band employing the on-the-fly (OTF) mode and producing maps of sizes typically of 1′× 1′. Using double-sideband receivers, the IRAM 30 m data cover a broader range of frequencies than the NOEMA data, between ≈ 213 and ≈ 221 GHz in the lower sideband and between ≈229 and ≈236 GHz in the upper sideband. At 1.37 mm, the angular and spectral resolution of the single-dish observations is ≈ 11′′ and 0.195 MHz (or 0.27 km s⁻¹), respectively. Merging the IRAM 30 m telescope data and the NOEMA visibilities yields coverage of the uv-plane in the inner 15 m. The IRAM 30 m data were calibrated using the CLASS package in gildas. To match the NOEMA observations, the IRAM 30 m data were smoothed to a spectral resolution of 3 km s⁻¹. A more detailed description of the calibration of the single-dish data is provided in Mottram et al. (2020).

2.2 POETS

The POETS survey complements multi-epoch VLBI observations of water masers with high-angular resolution deep imaging of radio continuum emission with the JVLA in a relatively large sample (38) of massive protostellar outflows. Moscadelli et al. (2016) give a general description of the POETS program, sample selection, observations, and data analysis. Hereafter, we resume only the main observational parameters of IRAS 21078+5211.

2.2.1 VLBI of 22 GHz water masers

The water maser emission (at 22.23508 GHz) in IRAS 21078+5211 was observed within the VLBA BeSSeL⁵ survey (project code: BR145E) at six epochs spanning about one year (from May 2010 to May 2011) with a totalobserving time per epoch of approximately seven hours. A general description of the BeSSeL VLBA observational setup is given in Reid et al. (2009). In the course of the observations of IRAS 21078+5211, the FWHM size of the VLBA beam varied from epoch to epoch in the range 0.5–0.9 mas, and the instantaneous field of view, limited by time-smearing, was ≈ 5.′′4. The velocity resolution for the 22 GHz water maser data was 0.42 km s⁻¹. In IRAS 21078+5211, maser absolute positions and proper motions are derived with an accuracy of 2 mas and 1–3 km s⁻¹, respectively.

2.2.2 JVLA continuum observations

IRAS 21078+5211 was observed using the JVLA of the National Radio Astronomy Observatory (NRAO⁶) between October 2012 and January 2013 in A configuration (project code: 12B-044), in C, Ku, and K bands (centered at 6.2, 13.1, and 21.7 GHz, respectively). In the following, we also refer to C- and K-band observations as the JVLA 5 and 1.3 cm continuum, respectively. The total on-source time was 15 min at C and Ku bands, and 30 min at K band. We employed the capabilities of the new WIDAR correlator to record dual polarization across a total bandwidth per polarization of 2 GHz. For the JVLA A-Array continuum images of IRAS 21078+5211, the FWHM size of the round synthesized beam and the largest recoverable angular scale are 0.′′ 44 and 4.′′ 5, 0.′′ 17 and 1.′′ 8, 0.′′ 08 and 1.′′ 2 at C, Ku and K band, respectively. The absolute positions of the images are expected to be accurate to within 30–40, 20, and 10 mas at 6, 13, and 22 GHz, respectively. The rms noise level of the continuum images is 10.7, 9.2, and 13.9 μJy beam⁻¹, at C, Ku, and K band, respectively.

2.3 LBT

Near-infrared (NIR) observations toward IRAS 21078+5211 were carried out with the instrument LUCI 2 (Ageorges et al. 2010; Seifert et al. 2010) at the 2 × 8.4 m LBT (Mount Graham, Arizona) on June 12, 2017, as part of the “Istituto Nazionale di Astrofisica” (INAF) program 2016_2017_54 (PI: Moscadelli). A field roughly centered on the position of IRAS 21078+5211 was imaged through the broad-band K_s (λ_c = 2.16 μm, FWHM = 0.27 μm) and narrow-band H₂ (λ_c = 2.127 μm, FWHM = 0.023 μm) filters. We operated LUCI 2 in seeing-limited imaging, single telescope mode, using camera N3.75, which provides a plate scale of ≈0.12′′ per pixel and a field of view of 4′ × 4′. The images were taken according to a dithering pattern alternating a pointing in which the nominal position of IRAS 21078+5211 is only a few arcseconds from the frame center and a pointing in which this position is ≈1′ from the frame center, so that all the frames contain the target position and skies can be obtained by selecting subsets of the more distant pointings. The H₂ observations consist of 31 dithered images with DIT = 6 s, NDIT = 25 (total integration time 4650 s) and the K_s observations consist of 13 dithered images with DIT = 2.74 s, NDIT = 18 (total integration time 641 s).

Data reduction was performed using standard IRAF⁷ routines. Each frame was flat-fielded using dome flat images and corrected for bad pixels. As they are vignetted, we cut out the outermost pixels before further steps. Skies were then constructed for each frame by median-filtering a subset of four frames selected from the nearest in time that also have large pointing offsets. For each filter, all sky-subtracted frames were registered and averaged together. The final mosaicked images exhibit point spread functions (PSFs) ≈1′′ wide. The continuum emission in the H₂ filter was estimated by comparing photometry of stars in the K_s and H₂ images. The two images were then scaled and subtracted from each other to produce an image showing pure H₂ 2.12 μm line emission.

As the images showed a complex blob of line emission near the target (which in turn appears obscured in the NIR), we decided to perform further adaptive optics (AO) assisted observations to try resolving possible bow-shock features making up the blob. New NIR observations were then obtained with LUCI 1 and FLAO (Esposito et al. 2012) at the LBT in diffraction-limited mode on October 1, 2017, as part of the INAF program 2017_2018_36 (PI: Massi). A field including both the target nominal position and the line emission blob was imaged through the H₂ and K_s filters using camera N30, which provides a plate scale of ≈ 0.015′′ per pixel and a field of view of 30′′ × 30′′. A star with R ≈ 12 at a distance of ≈1′ from the imaged field was selected as a guide star. A dithering pattern with random pointings a few arcseconds apart was selected. The H₂ data consist of a set of 24 dithered frames with DIT = 10 s and NDIT = 15 (total integration time 1 h) and the K_s data of a set of 6 dithered frames with DIT = 5 s and NDIT = 20 (total integration time 10 min). The frames were flat-fielded and bad pixels were removed. For each frame, a sky was constructed by median filtering the nearest (in time) four frames. The sky-subtracted frames were finally registered and averaged together. Unfortunately the AO correction was not very efficient (owing to the lack of a suitable guide star) and a Strehl ratio of only ≈ 1% was obtained. Nevertheless, these latest images (with a PSF FWHM ≈ 0.′′3) exhibit a better spatial resolution than the previous seeing-limited images. An image showing pure H₂ 2.12 μm line emission was obtained following the same procedure as above.

Comparing NIR images with radio and millimeter interferometric data, which usually have high positional accuracy, requires an astrometric calibration as precise as possible. First, we obtained an astrometric solution for the seeing-limited, large field-of-viewimages by cross-matching them with 2Mass Point Source Catalog entries. We were able to use about 100 relatively bright stars all over the field. Then, we cross-correlated the seeing-limited images and the AO-assisted images, finding ≈20 suitable stars to obtain an astrometric solution for the latter ones. We estimate that the astrometric accuracy of the AO-assisted images is better than 0.′′1.

3 Results

3.1 The cluster

The upper panel of Fig. 1 shows the 1.37 mm continuum emission in IRAS 21078+5211, which is already presented in Beuther et al. (2018, see Fig. B.3, upper panels). Using the CLUMPFIND algorithm, Beuther et al. (2018, see Table 5) identify 20 compact cores with intensities in the range 3–35 mJy beam⁻¹, the seven strongest of which are indicated in Fig. 1 (upper panel). The global pattern of the 1.37 mm continuum is elongated in the SE-NW direction: the PA of the major axis, determined through a linear fit to the positions of the seven strongest cores, is 149° ± 17°. The three strongest cores are found inside a higher-emission plateau at the center of the cluster, and the lower-intensity cores draw an arc of weaker emission encircling the plateau clockwise from SE to NW, with no cores found to SW.

The spatial distribution of the cores in the cluster revealed by the 1.37 mm continuum emission on linear scales of 10⁴ au has several properties in common with the patterns of different SF indicators on larger scales. The lower panel of Fig. 1 shows that the core cluster is embedded within a molecular cloud traced by the 850 μm emission from cold dust, which presents a similar SE-NW elongation (PA ≈ 140°). In turn, this molecular cloud is edged on the northeastern side with a necklace of compact 4.6 μm sources, associated with less embedded and, likely, more evolved YSOs in the region. This distribution resembles that of the cores in the cluster on ten times smaller scales.

To study the kinematics of the gas in the parental molecular cloud enshrouding the core cluster, we use the IRAM 30 m maps of three low-excitation (E_u ∕k_B ≤ 21 K) molecular lines: ¹³CO J = 2–1, C¹⁸O J = 2–1, and H₂ CO $J_{K_a,K_c}=$3_0,3–2_0,2. Figure 2 shows that, close to V_sys = −6.1 km s⁻¹, the velocities increase regularly from NW to SE along the major axis of the molecular cloud, reaching the highest values in correspondence of the core cluster. To investigate the gas kinematics inside the cluster, Fig. 3 shows velocity-channel maps of the merged NOEMA–IRAM 30 m observations of the ¹³CO J = 2–1 transition. In the two channels closer to V_sys, namely at V_LSR of −6.5 and −3.5 km s⁻¹, the ¹³CO J = 2–1 emission is prominent and extended on scales of a few 10⁴ au along a SE-NW direction close to the similarly oriented, elongation axes of the 1.37 mm core cluster and 850 μm molecular cloud. Despite the limited velocity resolution of only 3 km s⁻¹, a V_LSR gradient across the whole cluster, and its gaseous envelope, is clearly detected, the gas to N-NW emitting at lower V_LSR than the gas to SE. Thus, the merged NOEMA–IRAM 30 m observations of the ¹³CO J = 2–1 line confirm the velocity gradient observed at larger scales using only the IRAM 30 m data.

At high velocities, the spatial distribution of the more diffuse gas inside the cluster is revealed by the ¹³CO J = 2–1 (Fig. 3) and SO J_N = 6₅–5₄ (Fig. 4) transitions. At relatively high absolute line of sight (LOS) velocities, that is |V_LSR − V_sys|≥ 3 km s⁻¹, the channel emission in both lines is characterized by a SW-NE elongated feature pointing to core 1 (the strongest 1.37 mm continuum emitter), and placed either NE or SW of core 1 at high red-shifted (≥ −2 km s⁻¹) or blue-shifted (≤ −10 km s⁻¹) V_LSR. A straightforward interpretation for this emission feature in the ¹³CO and SO lines, as discussed in detail in Sect. 4.2.1, is in terms of a collimated outflow emitted by a YSO inside core 1. In the SO J_N = 6₅–5₄ channel maps at V_LSR ≈ −13 km s⁻¹, the emission is dominated by a strong bow shock of this outflow located ≈ 4′′ SW of core 1. Aside from the prominent signature of the outflow from core 1, the only other notable emission feature at high absolute LOS velocities is a compact source visible (especially in the SO channel maps) at very blue-shifted velocities (−28 km s⁻¹ ≤ V_LSR ≤ −20 km s⁻¹) to N-NE of core 4 (the northernmost core).

In IRAS 21078+5211, the higher-excitation (E_u∕k_B ≥ 70 K) molecular lines observed in the CORE program are only detected from inside the higher-emission plateau at the center of the 1.37 mm continuum distribution. The upper panel of Fig. 5 shows the velocity-integrated intensity of the CH₃CN J_K = 12₃–11₃ line, which is representative of the spatial distribution of all the observed dense-gas tracers. The CH₃CN emission peaks at the position of core 1 and it becomes progressively weaker moving from core 3 (near core 1) to core 2 (displaced further to S; see Fig. 1, upper panel). In the lower panel of Fig. 5, the 1.37 mm continuum is overlaid with the tracers of the outflows discovered in this region, that is the SW-NE (PA = 43° ± 10°), double-component radio jet detected in the POETS survey and the red- and blue-shifted integrated emission of the CORE SO J_N = 6₅–5₄ line. It should be noted that the NE component of the radio jet with thermal emission from ionized gas near the YSO, coincides with core 1, and that the radio jet and the SO outflow are approximately parallel to each other. These two findings lead us to assume that the YSO embedded in core 1 is driving the radio jet, which, in turn, powers the larger-scale molecular outflow. In Sect. 4.2.1, we check this assumption by comparing the properties of the radio jet and the SO outflow. Hereafter, we refer to the YSO inside core 1 as “YSO-1”.

Fig. 1 Upper panel: NOEMA 1.37 mm continuum emission (color map) of IRAS 21078+5211. The white contours correspond to levels ranging from 3σ (σ = 0.5 mJy beam−1) to the peak intensity of 35 mJy beam−1, in steps of 7σ. The little black crosses indicate the positions of the seven strongest 1.37 mm peaks labeled as in Beuther et al. (2018, see Table 5), with label numbers ordered in intensity. The restoring beam of the 1.37 mm continuum map is shown in the lower left corner. Lower panel: Wide-field Infrared Survey Explorer (WISE) 4.6 μm image (color map, Wright et al. 2010) toward IRAS 21078+5211. The white contours represent the SCUBA 850 μm continuum (Di Francesco et al. 2008), showing 10 levels increasing logarithmically from 7σ (σ = 0.1 Jy beam−1) to 11.4 Jy beam−1. The magenta rectangle delimits the region plotted in the upper panel.

Fig. 2 IRAM 30 m data. Intensity-weighted velocity (color map) of the C18O J = 2–1 (upper panel), 13CO J = 2–1 (middle panel), and H2CO $J_{K_a,K_c}=$30,3–20,2 (lower panel) lines determined over a small velocity range, [−7.3, −3.5] km s−1, around Vsys = −6.1 km s−1. The plotted region corresponds to the area in which the velocity-integrated emission of the C18O J = 2–1 line is higherthan 44 Jy beam−1 km s−1, ≈10σ. In each panel, the black contours represent the SCUBA 850 μm continuum (Di Francesco et al. 2008), showing 10 levels increasing logarithmically from 1 to 11.4 Jy beam−1. The white rectangle delimits the region plotted in the upper panel of Fig. 1.

Fig. 3 Merged NOEMA–IRAM 30 m data. Each panel presents the emission of the 13CO J = 2–1 line (color map) in a different velocity channel. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted intensity (in square parentheses in milliJansky per beam) are reported. The red boxes and ticks identify the panels corresponding to the central velocities, that is |VLSR − Vsys| < 3 km s−1. The white crosses mark the positions of the seven strongest peaks of the 1.37 mm continuum emission. At the central velocities, the black dashed line gives the major axis of the 1.37 mm core cluster.

3.2 Cluster center: cores 1, 2, and 3

Toward cores 1, 2, and 3, the emission of the higher-excitation molecular lines is strong enough to allow us to study the kinematics and physical conditions of the gas. For this purpose, we used the XCLASS (eXtended CASA Line Analysis Software Suite) tool (Möller et al. 2017). This tool models the data by solving the radiative transfer equation for an isothermal homogeneous object in local thermodynamic equilibrium (LTE) in one dimension. It can produce maps of column density, temperature, velocity, and line width by extracting spectra pixel-by-pixel from the image cubes and fitting all the unblended lines of a given molecular species simultaneously. The finite source size, dust attenuation, and line opacity are considered as well, which also permitsus to properly fit the transitions of relatively lower-excitation energy and higher optical depths (as the CH₃CN J_K = 12_K–11_K, K = 0–2 lines).

To complement the XCLASS results, we also analyzed several, higher-excitation molecular lines, individually.We verified that the selected lines are sufficiently optically thin to reliably trace the more embedded gas kinematics. More specifically, in this analysis we made use of the following molecular lines: CH₃CN J_K = 12_K–11_K (K = 3–6), CH₂CO $J_{K_a,K_c}=$11_1,11–10_1,10, and HC₃N J = 24–23. We fitted the line profile of these transitions toward the 1.37 mm continuum peak in core 1 and found that the optical depth is always < 0.5 at the peak velocity and decreases to values ≲0.05 across the line wings. The good quality of the fit for the CH₃CN J_K = 12_K-11_K (K = 0–6) lines is illustrated in Fig. 6.

We used XCLASS to simultaneously fit the emission of the CH₃CN J_K = 12_K–11_K (K = 0–6) lines (including the CH₃ ¹³CN isotopologs). Figure 7 shows the maps of the CH₃CN column density and velocity toward the central plateau of the 1.37 mm continuum, where the three most intense cores reside. The CH₃CN column density reaches 10¹⁶ cm⁻² toward core 1 and it is a factor of 2 lower toward core 2. Two V_LSR gradients are clearly detected: the first across the adjacent cores 1 and 3, characterized by a change of ≈ 4 km s⁻¹ over ≈ 1500 au; the second across core 2, with a change of ≈ 3 km s⁻¹ over ≈ 1000 au. In the following, we determine the velocities of the cores 1 and 3 to assess the contribution of their relative motion to the observed velocity gradient.

While the emission of the CH₃CN J_K = 12_K–11_K (K = 2–6) lines (with E_u∕k_B varying in the range 97–326 K) is concentrated mainly toward core 1 (see Fig. 7, upper panel), we also identified a few molecular lines of slightly lower excitation, such as the CH₂CO $J_{K_a,K_c}=$11_1,11–10_1,10 transition (E_u∕k_B = 76 K), which have comparable intensity toward cores 1 and 3 and can be employed to resolve the emission of the two cores in velocity. The upper panels of Fig. 8 show the maps of the peak emission of the CH₃CN J_K = 12₅–11₅ line toward core 1 and the CH₂CO $J_{K_a,K_c}=$11_1,11-10_1,10 line toward core 3, while the lower panel of Fig. 8 reports the spectra of the two molecular lines integrated over the corresponding core. Among the CH₃CN transitions, we selected the K = 5 component because it is intense and has a relatively high-excitation energy (E_u ∕k_B = 247 K), which allows us to better resolve the warm gas inside core 1. To determine the areas of the two cores, we used the minimum level (18.6 mJy beam⁻¹, indicated bythe black contours in Fig. 8) of the 1.37 mm continuum in correspondence of which two disconnected emission islands are still obtained. We stress that in this analysis we used the emissions of the CH₃CN J_K = 12₅–11₅ and CH₂CO $J_{K_a,K_c}=$11_1,11–10_1,10 lines only to derive the velocities of the cores 1 and 3, whose existence and geometrical properties are independently established from the 1.37 mm continuum emission. The two cores are separated by ≈ 770 au and their LOS velocities differ by ≈ 3 km s⁻¹. Therefore, it appears that the relative motion of the two cores can account for the observed V_LSR gradient. This result is used in Sect. 4.1.2, where we consider the possibility that the two cores are members of a binary system.

Our recent H₂ 2.12 μm LBT observations toward IRAS 21078+5211 are shown in Fig. 9 (upper panel) together with the radio continuum and SO line emissions tracing the jet from YSO-1 and the molecular outflow from core 1, respectively. The H₂ 2.12 μm emission arises from shock-excited, compressed, and hot molecular gas placed ≈ 2′′ SW of core 1, pinpointed by the NE component of the radio jet at the center of the SO molecular outflow. Examining the lower panel of Fig. 9, which provides a higher-angular resolution view of the shock structure, we note two well-shaped bow shocks at the SW tips of the H₂ 2.12 μm emission. Since the major axes of these bow shocks approximately coincide with the directions to core 1, it is plausible that these shocks have been excited by an outflow emerging from core 1. Assuming that the outflow source resides in core 1, Fig. 9(upper panel) shows that different outflow tracers, namely, the nearby SW lobe of synchrotron emission, the weak 5 cm spur ≈ 7000 au to NE, and the H₂ 2.12 μm bow shock ≈40 000 au to SW, are oriented at different PAs. In Sect. 4.4, we discuss several scenarios to explain the spread in the directions of the various outflow tracers.

Fig. 4 Merged NOEMA–IRAM 30 m data. Each panel presents the emission of the SO JN = 65–54 line (color map) in a different velocity channel. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted intensity (in square parentheses in milliJansky per beam) are reported. The red boxes and ticks identify the panels corresponding to the central velocities, that is |VLSR − Vsys| < 3 km s−1. The white crosses mark the positions of the seven strongest peaks of the 1.37 mm continuum emission.

3.3 Most prominent core 1

Toward IRAS 21078+5211, the linear resolution achieved with the CORE data is ≈625 au (for a FWHM beam size of ≈0.′′38). As seen in Sect. 3.2, this is sufficient to resolve nearby molecular cores, on scales ≳1000 au. Star formation models predict that the gravitational pull of high-mass YSOs influences the gas kinematics mainly at distances ≤1000 au (see, for instance, Kölligan & Kuiper 2018), where kinematic structures such as accretion disks and disk winds (DW) are expected. Toward core 1, the high signal-to-noise ratio of the CH₃CN J_K = 12_K–11_K (K = 2–6) and other high-excitation molecular lines (such as OCS J = 18–17, HC₃N J = 24–23, and HNCO $J_{K_a,K_c}=$10_1,10–9_1,9) permits us to study the gas kinematics on linear scales ≤1000 au. For each of these strong molecular lines, we fit a 2D Gaussian profile to the emission in each velocity channel and determine the peak position as a function of velocity. If the emission is unresolved, the positional accuracy is equal to (δθ ∕ 2) × (σ ∕ I) (see, e.g., Reid et al. 1988), where δθ is the FWHM beam size and I and σ are the peak intensity and the rms noise, respectively, of a given channel. With an average value of the ratio (I ∕ σ) = 27, the error is ≲10 mas. For channels with resolved emission, the centroid positions convey more complex kinematic information because in these channels it is very likely that we are observing the combination of different types of motions and the shape of the emission can strongly deviate from a Gaussian. Consequently, for these channels the peak position obtained with the Gaussian fit is less reliable.

The upper panel of Fig. 10 shows the distribution of channel centroids for the CH₃CN J_K = 12_K–11_K (K = 3–6) and HC₃N J = 24–23 lines, which are the most intense and compact molecular emissions in core 1. In Fig. A.1, the distribution of channel centroids is plotted separately for the CH₃CN and HC₃N lines to show that the two molecules have a very similar V_LSR pattern. This allows us to combine their emissions in our kinematical analysis. The colored contours of Fig. 10, reproducing the half-peak emission levels, clearly indicate that the structure is really compact only at the most extreme red- and blue-shifted velocities. At V_LSR ≈ −4 km s⁻¹, the contamination from the nearby core 3 is evident. In the lower panel of Fig. 10 only the positions of the compact-emission channels at the extreme velocities are plotted. The distribution is bipolar and elongated along a direction, at PA = 123° ± 23°, approximately perpendicular to the radio jet. The compact 1.3 cm continuum emission, which best pinpoints the position of YSO-1, is located at the center of the distribution. In Sect. 4.2, we propose that the derived velocity pattern traces the accretion disk around YSO-1.

Fig. 5 Upper panel: NOEMA data. Velocity-integrated intensity (color map) of the CH3 CN JK = 123–113 line. The 1.37 mm continuum emission is represented with black contours, showing seven levels increasing logarithmically from 11 to 35 mJy beam−1. The positions of the three strongest peaks of the 1.37 mm continuum emission are indicated, using the same labels as in Beuther et al. (2018, see Table 5). The restoring beams of the CH3 CN JK = 123–113 line and 1.37 mm continuum maps are shown in the lower right and left corners, respectively. Lower panel: the blue and red contours reproduce the emission of the SO JN = 65–54 line integrated over the velocity ranges [−22,−13] and [−1, 5] km s−1, respectively. The plotted levels are from 0.03 to 2.5 Jy beam−1, in steps of 0.27 Jy beam−1, and from 0.6 to 2.2 Jy beam−1, in steps of 0.37 Jy beam−1, for the blue- and red-shifted SO emission, respectively. The magenta contours give the JVLA A-Array continuum at 5 cm (Moscadelli et al. 2016), showing levels at 40, 50, and 90% of the peak emission of 95 μJy beam−1: the magenta dashed line connects the two strongest (nearby, but resolved) 5 cm peaks. The two black dashed lines delimit the viewing angle of the blue-shifted emission of the SO JN = 65–54 line from the NE 5 cm peak (aligned in position with core 1, see Sect. 3.1). The black contours have the same meaning as in the upper panel.

Fig. 6 Beam-averaged spectrum (black line) of the CH3CN JK = 12K–11K (K = 0–8) line emissiontoward the 1.37 mm continuum peak. The positions of the K = 0–8 lines are indicated. The fits (red histograms) of the K = 0–6 lines performed with XCLASS are also shown.

4 Discussion

4.1 The cluster

4.1.1 Mass flow toward the core cluster

The 1.37 mm continuum cluster is found at the density peak of a SE-NW elongated molecular cloud edged on the northeastern side with a necklace of less embedded YSOs and on the southwestern side with a filament of diffuse 4.6 μm emission (see Fig. 1, lower panel). These findings suggest that the molecular cloud is a density enhancement of a more extended, infrared dark filament. In Sect. 3.1 we showed that the velocities of the slow-moving (≤ 3 km s⁻¹) gas of the cloud and the cluster vary smoothly along the major axis of the molecular cloud over linear scales of 10⁴ –10⁵ au (see Figs. 2 and 3). This regular change in V_LSR with position makes us favor the interpretation in terms of a flow in the molecular gas (converging toward the density peak where the most massive cores reside) compared with cloud-cloud collision, which would rather appear as a sudden jump in velocity. In the assumption of a mass flow, the observed V_LSR gradient can be easily explained if the blue-shifted NW side of the molecular cloud is farther away from us than the red-shifted SE side. The spatial and velocity extents of this V_LSR gradient are comparable with those of the gradients, ~1 km s⁻¹ per 0.1 pc, observed in infrared dark clouds (IRDC) at the sites of YSOs (Ragan et al. 2012), which are often interpreted as infall. It is also consistent with simulations of large-scale accretion flows along filaments, gravitationally accelerated toward local density peaks (see, e.g., Tobin et al. 2012; Smith et al. 2013).

Assuming we are observing a mass infall toward the core cluster, we can use both the IRAM 30 m and merged NOEMA–IRAM 30 m observations of the ¹³CO J = 2–1 line to calculate the mass infall rate. The latter can be derived from the ratio of the momentum, P_inf, to the length, L_inf, of the flow. In Appendix B, we describe the method to calculate the momentum of a molecular outflow, by employing the ¹³CO J = 2–1 emission. The same equations hold to derive the momentum of a molecular infall as well. We employ the same parameters used for the ¹³CO J = 2–1 molecular outflow (see Table B.1) with the following exceptions. First, the velocities to be considered are those within a small interval about the systemic velocity, that is, [−7, −5.5] km s⁻¹ for the IRAM 30 m observations (see Fig. 2), and [−6.5, −3.5] km s⁻¹ for the merged NOEMA–IRAM 30 m data (see Fig. 3). Second, the FWHM size of the map beam is B_max = B_min = 11.′′8, for the IRAM 30 m observations. It is remarkable that, using either the IRAM 30 m or the merged NOEMA–IRAM 30 m observations, we obtain very consistent values for the infall momentum P_inf cos(i_inf) ≈ 28 M_⊙ km s⁻¹, where i_inf is the inclination angle of the flow with respect to the LOS. Taking L_inf sin(i_inf) ≈ 6 × 10⁴ au, corresponding to the average value between the (sky-projected) flow lengths measured in the IRAM 30 m (8 × 10⁴ au, see Fig. 2) and NOEMA-IRAM 30 m (4 × 10⁴ au, see Fig. 3) maps, we finally derive a mass infall rate Ṁ_inf cot(i_inf) = P_inf cos(i_inf) ∕ L_inf sin(i_inf) ≈ 28 M_⊙ km s⁻¹ / 6 × 10⁴ au ≈ 10⁻⁴ M_⊙ yr⁻¹. Following the discussion in Appendix B, owing to the uncertainties in the excitation temperature of the ¹³CO J = 2–1 line and abundance ratio of ¹³CO, our estimateof the infall rate should be accurate to within a factor of 5. The derived value of infall rate falls at the lower end of the range oflongitudinal flow rates measured in filaments associated with high-mass star-forming regions (Chen et al. 2019; Treviño-Morales et al. 2019).

Fig. 7 NOEMA data. The color map reproduces the column density (upper panel) and the velocity (lower panel) of the CH3 CN emission determined by fitting with XCLASS the CH3CN JK = 12K–11K (K = 0–6) lines (including the CH3 13CN isotopologs), simultaneously. The plotted regions correspond to the areas in which the 1.37 mm continuum emission is higher than 13 and 16 mJy beam−1 for the column density and velocity plots, respectively. The 1.37 mm continuum emission is represented with black contours, showing six levels increasing logarithmically from 13 to 35 mJy beam−1: the positions of the three strongest 1.37 mm peaks are denoted using the same labels as in Beuther et al. (2018, see Table 5). The restoring beams of the 1.37 mm continuum and CH3 CN maps are shown in the lower left corner of the upper panel and lower right corner of the lower panel, respectively.

4.1.2 Temperature and mass distribution over the cluster center

Figure 11 shows the rotational temperature determined with XCLASS by simultaneously fitting the emission of the CH₃CN J_K = 12_K–11_K (K = 0–6) lines (including the CH₃ ¹³CN isotopologs). While toward core 1 the gas temperature is always ≥115 K and may reach 200 K, toward core 3 it varies in the range 90–130 K. In correspondence of core 2 the temperature map is less homogeneous, probably owing to insufficient signal-to-noise ratio, and has an average temperature that is intermediate between that of core 1 and core 3. By combining the maps of the 1.37 mm dust continuum and gas temperature, assuming that the dust emission is optically thin, we can make a reasonable estimate of the core masses. For consistency with Beuther et al. (2018), we adopt a dust absorption coefficient of 0.9 cm² g⁻¹ (Ossenkopf & Henning 1994) and a gas-to-dust mass ratio of 150 (Draine 2011). Over the areas of core 2 and the combined cores 1 and 3 (see Fig. 11), the mass is found to be 0.46 and 0.92 M_⊙, respectively.Inside the individual cores 1 and 3, whose corresponding areas are delimited by the black contours in the upper panels of Fig. 8, we obtain a gas mass of 0.42 M_⊙ and 0.12 M_⊙, respectively. These values represent the mass inside the cores within radii ≤500 au, and they are lower than the total core masses determined by Beuther et al. (2018, see Table 5) over radii ≥1000 au with a mean gas temperature of 66 K.

In Sect. 3.2, using the CH₃CN J_K = 12_K–11_K (K = 0–6) lines (see Fig. 7), we traced the kinematics of the gas inside the three cores in the cluster center. Now, we use this result to derive the dynamical masses of the cores and estimate the masses of the embedded YSOs. In Sect. 4.2.2 we determine the mass of YSO-1 (the YSO inside core 1) to be 5.6 ± 2.0 M_⊙. In Sect. 3.2, we showed that the V_LSR gradient across the adjacent cores 1 and 3 can be well explained in terms of the relative motion of the cores. Assuming that we are observing rotation seen almost edge-on, in gravito-centrifugal equilibrium, from the spatial separation between the cores, ≈770 au, and their difference in velocity, ≈3 km s⁻¹ (see Sect. 3.2), we derive a total mass of ≈8 M_⊙. Subtractingboth the mass of YSO-1 of ≈5.6 M_⊙ and that (0.92 M_⊙, see above) of the gas embedded inside cores 1 and 3, we obtain a rough estimate for the mass of the YSO inside core 3 of ≈1–2 M_⊙. Regarding core 2, interpreting the V_LSR gradient (≈3 km s⁻¹ across ≈1000 au, see Sect. 3.2) observed toward this core as rotation, in gravito-centrifugal equilibrium, we derive a dynamical mass of ≈2.5 M_⊙. After subtracting the gas mass of 0.46 M_⊙ (see above), the derived mass for the YSO inside core 2 is ≈2 M_⊙. Despite the large uncertainty in these determinations, the presence of relatively low-mass YSOs inside cores 2 and 3 agrees with the lower temperature (see Fig. 11) and mass of these cores with respect to core 1.

Fig. 8 NOEMA data. Upper panels: the color maps reproduce the emission of the CH3 CN JK = 125–115 (left) and CH2CO $J_{K_a,K_c}=$111,11–101,10 (right) lines. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted values (in milliJansky per beam) are reported. The white contours reproduce the 1.37 mm continuum, with the same levels as in Fig. 7: the positions of the two nearby 1.37 mm peaks are denoted using the same labels as in Beuther et al. (2018, see Table 5), namely, 1 and 3 for the primary and secondary peak, respectively. Right panel: the blue segment connects the two peaks. In the two panels, the black contours delimit the integration areas used to produce the spectra toward cores 1 and 3, presented in the lower panel. Lower panel: spectra of the CH3 CN JK = 125–115 line toward core 1 (blue histogram) and the CH2CO $J_{K_a,K_c}=$111,11–101,10 line toward core 3 (red histogram). The flux scales for the CH3CN and CH2CO spectra are reported on the left and right, vertical axes, respectively. The blue and red dashed vertical lines indicate the approximateVLSR of the cores 1 and 3, respectively.

4.1.3 Present state and future evolution

Among the 20 cores identified in the 1.37 mm continuum emission of IRAS 21078+5211 by Beuther et al. (2018, see Table 5), only cores 1, 2, and 3, which are the most massive (≥2 M_⊙), show signatures of an embedded protostar. As described in Sect. 3.1 (see also Fig. 5, upper panel), the higher-excitation molecular lines are only detected toward these three cores, indicating that they are sufficiently warm to require local heating by a protostar. Fitting the CH₃CN J_K = 12_K–11_K lines, the gas temperature over the three cores is found everywhere ≳100 K (see Sect. 4.1.2), and, through the observed V_LSR gradients (see Sect. 3.2), the masses of the embedded YSOs are evaluated in the range 1–6 M_⊙ (see Sect. 4.1.2). The concentration of the most massive cores (and YSOs) toward the center of the 1.37 mm continuum cluster hints at mass segregation.

Judging the evolutionary state of the other, less massive cores is difficult with the present observations. Comparing with the most massive cores, the non-detection toward these cores of the higher-excitation molecular lines indicates an upper limit of ≲ 1 M_⊙ for the mass of the embedded YSOs. The fact that no molecular outflows are detected from the less massive cores (see Sect. 3.1) would be consistent with most of these cores still being in a prestellar phase. However, with our data, a molecular outflow is observed only from the most massive YSO, YSO-1 in core 1, and this hints at a sensitivity limit when using the ¹³CO J = 2–1 and SO J_N = 6₅–5₄ lines as outflow tracers. Future interferometric observations in the CO J = 2–1 line, more intense and with broader excitation conditions, could unveil less massive and energetic outflows. On a more speculative ground, the non-detection of outflows from the less massive cores could be also explained if the (putative) embedded low-mass YSOs were not actively accreting because most of the gas available across the cluster was gathered by YSO-1.

To discuss the future evolution of the IRAS 21078+5211 core cluster is useful to compare it with more evolved states as the well-studied, NIR stellar clusters around Ae or Be pre-main-sequence (PMS) intermediate-mass (1–6 M_⊙) stars (see, for instance, Testi et al. 1999). Ae and Be PMS stars are young enough, 0.5–5 Myr, that any population of lower-mass stars born in the same environment did not have enough time to move away from the birthplace. Therefore, the spatial distribution of the stars reflects that of the parental molecular cores. The typical size of these NIR stellar clusters is ≈0.2 pc, and the highest stellar densities, up to 10³ pc⁻³, are found in correspondence with the most massive (≈6 M_⊙) Ae or Be stars. In IRAS 21078+5211, 20 cores are detected over an area of ≈ 10⁻² pc², equivalent to a density of 2 × 10⁴ pc⁻³, and the most massive YSO (YSO-1) has already attained a mass close to 6 M_⊙. In Sect. 4.1.1, we calculated a mass infall rate of ~10⁻⁴ M_⊙ yr⁻¹ from the parental molecular cloud (total mass of 177 M_⊙, Beuther et al. 2018, see Table 1) to the core cluster, which, over the characteristic formation time of a few 10⁵ yr of the Ae and Be PMS stars (Palla & Stahler 1993), would imply an infall toward the cluster of a few 10 M_⊙. If a relevant fraction of this mass is accreted by the most massive YSOs at the center of the cluster, it is very likely that at least one high-mass (≥8 M_⊙) star will form. In conclusion, because of the relatively high stellar density and the likelihood of ultimately forming massive stars, we think that the IRAS 21078+5211 core cluster can be the progenitor of a stellar system whose most massive members are early B or late O-type high-mass stars, rather than Ae or Be intermediate-mass stars.

Fig. 9 Comparison of CORE, POETS and LBT observations. Upper panel: the gray-scale map reproduces the H2 2.12 μm emission observed with LBT toward IRAS 21078+5211. Referring to Fig. 5, lower panel, the blue and red contours have the same meaning, and the yellow contours replace the magenta contours for the JVLA 5 cm continuum. The magenta dashed, white dotted, and cyan dashed lines connect the main 5 cm peak with the secondary (nearby) 5 cm peak, the weak 5 cm spur to NE, and the bow shock traced by the LBT H2 2.12 μm emissionto SW, respectively. The white dashed box delimits the region plotted in the lower panel. Lower panel: adaptive-optics assisted LBT observations of the H2 2.12 μm emission (gray-scale map): the cyan dashed line has the same meaning as in the upper panel.

Fig. 10 NOEMA data. Upper panel: the gray-scale map reproduces the 1.37 mm continuum emission, plotting values in the range 10–35 mJy beam−1. The colored dots give the (Gaussian-fitted) positions of the channel emission peaks for the CH3 CN JK = 12K–11K (K = 3–6) and HC3N J = 24–23 lines; the colors denote the channel VLSR. The colored contours indicate the half-peak level for the CH3CN JK = 123-113 emission in individual channels. The yellow contours show the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels at 70, 80, and 90% of the peak emission of 0.50 mJy beam−1. Lower panel: the gray-scale map, yellow contours, and colored dots have the same meaning as in the upper panel. Only the channel emission peaks at the extreme blue- and red-shifted velocities (reported in the lower right corner of the panel) are shown. The white continuous and dashed lines indicate the best-fit major axis of the distribution of the peaks and the corresponding fit uncertainty, respectively.

Fig. 11 NOEMA data. The color map shows the distribution of the rotational temperature derived by fitting with XCLASS the emission of the CH3CN JK = 12K–11K (K = 0–6) lines (including the CH3 13CN isotopologs). The plotted region corresponds to the area in which the 1.37 mm continuum emission is higher than 16 mJy beam−1. The black contours reproduce the 1.37 mm continuum, showing the same levels as in Fig. 7: the positions of the three strongest 1.37 mm peaks are indicated using the same labels as in Beuther et al. (2018, see Table 5). The restoring beam of the CH3 CN maps is shown in the lower right corner of the panel.

4.2 Disk-jet system associated with YSO-1

4.2.1 Radio jet and molecular outflow

Absolute positions and proper motions of the 22 GHz water masers observed in IRAS 21078+5211 with the POETS survey are presented in Fig. 12. The water masers show a SW-NE elongated distribution at radii between 100 au and 300 au from the peak of the compact 1.3 cm JVLA continuum, corresponding to the core of the radio jet and the YSO-1 position. The elongation, at PA ≈ 38°, of the maser spatial distribution, and the average direction of the proper motions, PA ≈ 49° (Moscadelli et al. 2019,see Table A1), are in good agreement with the orientation, PA ≈ 43°, of the double-component radio jet (see Fig. 5). The faint, thermal, and nonthermal emissions from the NE and SW radio lobes, respectively, and the water masers are all manifestations of different types of shocks produced by the jet. In the jet core, close to YSO-1, we observe free-free emission, which likely originates from relatively weak (C-type), internal shocks of the jet corresponding to changes in the mass loss or ejection velocity (Moscadelli et al. 2016; Anglada et al. 2018); the nonthermal emission is likely synchrotron emission from electrons accelerated at relativistic velocities in strong jet shocks through first-order Fermi acceleration (Bell 1978); finally, the fast water masers emerge from relatively strong (J-type) shocks, whereas the jet impinges on very dense (molecular hydrogen number density, $n_{{\text{H}}_2}~{10}^8$ cm⁻³) ambient material at high velocity (Moscadelli et al. 2020).

As reflected in Fig. 5 (lower panel), the radio jet and the molecular outflow traced by the SO J_N = 6₅–5₄ line at larger scales are approximately parallel to each other and, as already anticipated in Sect. 3.1, the molecular outflow could be powered by the jet. To assess the physical association between the jet and molecular outflow, in the following we determine and compare their physical properties. In Appendix B, we describe the method to calculate the momentum of the molecular outflow, P_out. We obtain P_out cos(i_out) ~ 50 M_⊙ km s⁻¹, where i_out is the inclination angle of the outflow with respect to the LOS. Employing the water maser positions and three-dimensional (3D) velocities, the momentum rate of the jet in IRAS 21078+5211 was estimated to be Ṗ_jet ~ 1 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹ by Moscadelli et al. (2016, see Sect. 6.2, Eq. (1), and Table 5). In this calculation, a major source of uncertainty was the pre-shock, ambient density, which was assumed to be $n_{{\text{H}}_2}~$ 10⁸ cm⁻³. Now, thanks to the CORE data, we can better constrain this parameter. Averaging the derived mass distribution (see Sect. 4.1.2) over a cylindric volume equal to the beam area times the beam size, in the direction of the continuum peak we find the value $n_{{\text{H}}_2}\approx$ 3 × 10⁸ cm⁻³, which should be accurate to within a factor of a few. Accordingly, our new, more precise, estimate of the momentum rate of the jet is Ṗ_jet ~ 3 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹.

We wish to derive the momentum rate of the molecular outflow, to be compared with that of the jet. Knowing the value of P_out cos(i_out) (see above), we need to estimate the inclination angle and timescale of the outflow. The clear spatial separation of the blue- and red-shifted lobes (see Fig. 5, lower panel) indicates that the gas flows along directions significantly inclined with respect to both the LOS and the plane of the sky. Based on the opening angle of ≈ 40° of the blue-shifted outflow lobe (see Fig. 5, lower panel), we judge that the outflow axis has to be at least 20° away from both the LOS and the plane of the sky, that is 20 ≤ i_out ≤ 70°. The dynamical timescale of the outflow, t_out, can be estimated from the ratio of the size of the lobes to the corresponding maximum V_LSR range. Referring to the SO J_N = 6₅–5₄ emission in Fig. 5 (lower panel), we have: t_out tan(i_out) ≈ 22′′ / 27 km s⁻¹ = 6.3 × 10³ yr, at a distance of 1.63 kpc. We caution that previous studies have pointed out that the dynamical timescale of molecular outflows can largely underestimate (by a factor between 5 and 10) the true outflow age (see, for instance, Parker et al. 1991). Accordingly, the upper limit for the momentum rate of the outflow is found to be Ṗ_out cos ²(i_out) / sin (i_out) ≲ 8 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹. Consideringthat this value is only accurate to within an order of magnitude and that 0.4 ≤ sin(i_out) / cos ²(i_out) ≤ 8 for 20° ≤ i_out ≤ 70°, we conclude that the momentum rate of the molecular outflow is consistent with that of the jet. It is then plausible that the jet powers the molecular outflow.

Fig. 12 POETS water maser observations. Upper panel: colored dots and arrows give absolute positions and proper motions of the 22 GHz water masers determined with VLBI observations (Moscadelli et al. 2016); colors denote the maser VLSR. The dot area scales logarithmically with the maser intensity. The black contours indicate the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels from 10 to 90% of the peak emission of 0.50 mJy beam−1 in steps of 10%. Lower panel: colored dots and black contours have the same meaning as in the upper panel. Colored cones are employed to visualize the maser 3D velocities, representing the inclination with respect to the LOS through the ellipticity of the cone basis and the uncertainty in the direction by means of the cone aperture.

4.2.2 Mass and luminosity of YSO-1

The velocitypattern in Fig. 10 (lower panel) is consistent with that expected for a rotating disk around a YSO. It is elongated approximately perpendicular to the radio jet, the blue- and red-shifted velocities lie on the two sides of the YSO, and its size of ≈ 400 au agrees with the model predictions for accretion disks around high-mass YSOs (see, for instance, Kölligan & Kuiper 2018). The V_LSR does not present a regular change with the position, but that likely results from the limited precision in tracing the velocity pattern using channel emission centroids and the insufficient angular resolution.

Making the assumption (to be checked a posteriori, see below) of Keplerian rotation, we determine the central mass, M_⋆, by minimizing the following χ² expression: $${\chi }^2={\sum }_j\frac{{[V_j-(V_{\star }\pm 0.74{(M_{\star }{\sin }^2(i_{\text{rot}}))}^{0.5}|S_j-S_{\star }{|}^{-0.5})]}^2}{{(\Delta V_j)}^2},$$(1)

where V_j and S_j are the channel V_LSR (in kilometer per second) and corresponding peak positions (in arcsecond), the index j runs over all the fitted channels, and the + and − symbols hold for red- and blue-shifted velocities, respectively. V_⋆ (in kilometer per second) and S_⋆ (in arcsecond) are the V_LSR and position of YSO-1, respectively. The central mass M_⋆ is given in solar masses. Indicating with i_rot the inclination of the disk rotation axis with respect to the LOS, the factor sin ²(i_rot) takes into account that the disk plane is seen at an angle 90 − i_rot from the LOS and the observed V_LSR corresponds to the rotation velocity multiplied by the factor sin(i_rot). To take into account the uncertainty on the velocity and that on the position, the global velocity error ΔV_j is obtained by summing in quadrature two errors: that on the velocity (taken equal to half of the channel width) and that obtained by converting the error on the offset into a velocity error through the function fitted to the data.

The V_LSR, V_⋆ = − 6.4 km s⁻¹, and position S_⋆ = 0.′′81 of YSO-1 are assumed equal to the average of the V_LSR and positions of all the channels. Thus, the term M_⋆ sin²(i_rot) is the only free parameter in Eq. (1). We allowed this term to vary over the range 2–14 M_⊙, which is consistent with the upper limit of ≈ 12 M_⊙ imposed by the bolometric luminosity of ≈10⁴ L_⊙ (Davies et al. 2011). Figure 13 reports the plot of the χ² vs. the free parameter; the blue dashed lines delimit the 1σ confidence interval following Lampton et al. (1976). This plot shows that we do find an absolute minimum of the χ² and the determined 1σ range for the values of the free parameter is M_⋆ sin²(i_rot) = 5.3 ± 1.6 M_⊙.

The fitted value of M_⋆ is much larger than the gas mass, 0.42 M_⊙, inside core 1, derived from the dust emission in Sect. 4.1.2. As noted in Sect. 4.1.2, this is the mass inside core 1 within a radius ≤500 au, which is comparable to the size of the disk structure traced by the CH₃CN and HC₃N lines around YSO-1 (see Fig. 10, lower panel). The YSO-1 envelope should extend at significantlylarger radii and be much more massive. The finding that the central mass is much greater than the mass of the disk supports our assumption of Keplerian rotation. The disk, traced by the bipolar V_LSR pattern of molecularemissions (see Fig. 10, lower panel), and the radio-maser jet are approximately perpendicular to each other in the plane of the sky. This suggests that the jet is directed close to the rotation axis of the disk, or i_jet ≈ i_rot (where i_jet is the inclination of the jet with respect to the LOS). From the angular distribution of the maser proper motions, reported in Fig. 12, the semi-opening angle of the jet is evaluated to be ≈ 18° by Moscadelli et al. (2016, see Table 5). Figure 12, lower panel, also shows that the water masers are both blue- and red-shifted, which indicates that the jet intersects the plane of the sky and the relation |i_jet− 90°|≲ 18° must hold. Following these considerations, we can write sin²(i_rot) ≈ sin²(i_jet) ≳ sin²(72°) = 0.9. Allowing for the variation of the factor sin²(i_rot) in the range 0.9–1, the 1σ interval for the fitted YSO-1 mass is M_⋆ = 5.6 ± 2.0 M_⊙.

By employing the relationship between bolometric luminosity, L_bol, and outflow momentum rate, Ṗ, for massive YSOs by Maud et al. (2015) (Log₁₀ [Ṗ ∕ M_⊙ km s⁻¹ yr⁻¹] = −4.8 + 0.61 Log₁₀ [L_bol ∕ L_⊙]), and our best estimate of the momentum rate Ṗ_jet ~ 3 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹ (see Sect. 4.2.1), we infer a luminosity of 5 × 10³ L_⊙ for YSO-1. This value is consistent with the luminosity of 1.3 × 10⁴ L_⊙ of the IRAS 21078+5211 region, which was derived by Molinari et al. (1996) using low-angular resolution IRAS data (~1′, Neugebauer et al. 1984) only⁸. Indeed, combining the far-infrared IRAS fluxes with the mid-infrared, Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010) and Midcourse Space Experiment (MSX; Egan et al. 2003) data at higher angular resolution (~10′′), Moscadelli et al. (2016, see Table 1) determine a lower luminosity of 5 × 10³ L_⊙, which agrees well with that inferred for YSO-1.

Fig. 13 χ2 values from the Keplerian fit (see Eq. (1)) vs. the free parameter M⋆ sin 2(irot). The horizontal black dashed line indicates the value ${\chi }_{\text{min}}^2$+1.0 = 4.12, which corresponds to the 1σ confidence level for the χ2-distribution with onefree parameter (Lampton et al. 1976). The two vertical blue dashed lines indicate the 1σ confidence interval for the best-fit value of the free parameter.

4.3 Magnetized and rotating jet

The radio-maser jet emerging from YSO-1 is unique in combining two properties that, to our knowledge, for the first time are observed associated with the same object: a lobe with nonthermal emission and jet rotation. We noted the first point several times through this article and this refers to the SW lobe of the double-component radio jet (see Fig. 5, lower panel), whose spectral index over the frequency range 6–22 GHz is ≤ −0.7 (Moscadelli et al. 2016, see Table 4 and Fig. 19). Observations of nonthermal components in jets from high-mass YSOs are still rare, although good examples were recently discovered thanks to high-angular resolution, sensitive radio survey (for instance, the POETS source G035.02+0.35, Sanna et al. 2019). The observational evidence for the second point, as already discussed in Moscadelli et al. (2019, see Sect. 3.2), is the V_LSR gradient of the water masers transversal to the jet direction. Figure 12 (lower panel) shows the monotonic change of the maser V_LSR along a SE-NW direction perpendicular to the SW-NE collimation axis of the proper motions, with red- and blue-shifted V_LSR toward SE and NW, respectively. The linear correlation between the maser LOS velocities, V_LOS = V_LSR − V_⋆, and positions projected perpendicular to the jet direction, R_per, is clearly illustrated in Fig. 14. A straightforward interpretation for the maser 3D velocity pattern is in terms of a composition of two motions: a flow along, and a rotation around, the jet axis.

The findings of nonthermal emission from a jet lobe and jet rotation are strong indications that magnetic fields play an important role in both launching and accelerating the jet. The interpretation of the nonthermal continuum in terms of synchrotron emission requires the presence of a magnetic field sufficiently ordered and strong to trap electrons in jet shocks, confine them even at relativistic velocities, and warrant the detection of their synchrotron signal. Jet rotation is considered a critical test for magneto-centrifugal (MC) DW (Blandford & Payne 1982; Pudritz et al. 2007), where the magneto-centrifugally launched jetextracts the excess angular momentum from the disk gas, which, then, can be accreted by the protostar. Comparing Figs. 10 and 12, we evince that the jet and the disk of YSO-1 rotate about axes approximately parallel to each other and have the same sense of rotation, that is, toward and away from us to NW and SE, respectively.This result supports our idea that the jet rotation stems from the disk through the magnetic leverage of a MC DW.

The measurement of the 3D velocities of the jet near the YSO through water maser VLBI observations permits a more quantitative comparison with the predictions of the MC DW theory. According to the MC DW models (Pudritz et al. 2007), the YSO wind is centrifugally launched along the magnetic field lines threading the accretion disk, it is accelerated to the Alfvén speed sliding along the rotating field lines, and then is magnetically collimated. A large body of different MC DW simulations have shown that, depending on the assumed magnetic field configuration, the final jet structure can be either collimated toward a cylinder or at wide angle (Pudritz et al. 2007). In the following, we compare our water maserobservations with wind models that recollimate within a cylindrical flow tube (see, for instance, the simulations in Staff et al. 2015 and Kölligan & Kuiper 2018). We also note that a model of a cylindrical rotating jet has been recently proposed by Burns et al. (2015) to interpret the 3D motion of the water masers observed with VERA (VLBI exploration of radio astrometry; Kobayashi et al. 2003) toward the massive YSO S235AB-MIR.

Let us indicate with r₀ the radius of the footpoint of the magnetic field line anchored to the disk and with r_A the distancefrom the rotation axis at which the wind velocity equals the Alfvén speed. At radii r_c ≥ r_A, the wind recollimates, and its velocity has two (dominant) components: the main component oriented along the jet axis, ϒ_axi, and the azimuthal component, ϒ_Φ, representing the rotation around the axis. Combining Eqs. (12) and (13) of Pudritz et al. (2007), the ratio of these two components can be related to the radius r_c of the cylindrical flow through the expression $${{\rm Y}}_{\text{axi}}/{{\rm Y}}_{\Phi }\approx \sqrt{2}{r}_{\text{c}}/r_{\text{A}}.$$(2)

The water masers originate in shocks that arise where the jet impinges on high-density material and propagate away along the jet direction at a speed, assuming momentum conservation, V = $\sqrt{{\rho }_{\text{jet}}/{\rho }_{\text{amb}}}$ ϒ (see, for instance, Masson & Chernin 1993). In this equation, ρ_jet and ϒ are the jet density and speed, respectively, and ρ_amb is the ambient density. Since the water masers move parallel to the wind, their motion can be described in terms of the composition of the same velocity components, that is, that directed along the jet axis and the azimuthal component. Let us indicate with V_axi and V_per the maser velocity components in the plane of the sky parallel and perpendicular to the sky-projected jet axis (at PA = 43°). Since the jet is directed close (≤18°, see Sect. 4.2.2) to the plane of the sky, V_axi approximates the component along the jet axis with a maximum error of cos ⁻¹(18°) ≈ 5%. The azimuthal component, V_Φ, is equally well approximated with $V_{\Phi }\approx \sqrt{V_{\text{per}}^2+V_{\text{LOS}}^2}$. From the considerations above, it is clear that the ratio of the maser velocity components V_axi and V_Φ is an accurate measurement of the corresponding ratio of the wind velocity components, that is, ϒ_axi ∕ ϒ_Φ ≈ V_axi ∕ V_Φ.

Figure 15 shows that the ratio V_axi ∕ V_Φ is approximately constant for the vast majority of the water masers (with measured proper motions). We obtain V_axi ∕ V_Φ ≈ 3.4 ± 1. The only notable exception is a small maser cluster at the most blue-shifted velocities (V_LOS ≤ −20 km s⁻¹) closer to YSO-1 (axis-projected separation R_axi ≤ 150 au), for which 1 ≲ V_axi ∕ V_Φ ≲ 2. On the basis of Eq. (2), a simple explanation for the small spread in the ratio of the maser velocity components is that most of the masers originate within a thin cylindrical shell, tracing shocked wind material MC-launched within a thin annulus of radius r₀ of the rotatingdisk. This interpretation fits with the models predicting a shock origin for the water masers, since it is plausible that the most external layer of the jet mainly interacts with the surrounding ambient material and the water masers are mainly found on the jet wall. Thus, a reasonable estimate for the radius of the jet, r_jet, can be obtained from half the spread in R_per, the maser positions projected perpendicular to the jet direction, considering only the water masers with V_LOS ≥ −20 km s⁻¹ (see Fig. 14). We derive r_jet ≈ 10 mas, or ≈16 au.

The terminal poloidal velocity of a MC DW, which for a cylindrical flow corresponds to ϒ_axi, can be approximated as (Pudritz et al. 2007, see Eq. (12)) $${{\rm Y}}_{\text{axi}}\approx \sqrt{2}\frac{r_{\text{A}}}{r_0}{\upsilon }_{\text{K},\text{0}},$$(3)

where υ_K,0 is the disk Keplerian velocity at the radius r₀. From Eq. (2) and knowing both the ratio of ϒ_axi / ϒ_Φ ≈ 3.4 and the radius r_c = r_jet ≈ 16 au, we get r_A ≈ 6.7 au. Numerical and theoretical models (for a review see Pudritz et al. 2007) supported by observations (Bacciotti et al. 2002) constrain the ratio of r_A ∕r₀ ≈ 3. Then, using the above value of r_A we infer r₀ ≈ 2.2 au. From the latter and the estimated mass of YSO-1 M_⋆ ≈ 5.6 M_⊙, we find υ_K,0 ≈ 50 km s⁻¹, which, through Eq. (3), yields ϒ_axi ≈ 200 km s⁻¹ and ϒ_Φ ≈ 200 km s⁻¹ / 3.4 ≈ 60 km s⁻¹.

The derived values of jet radius and terminal velocity are in reasonable agreement with the few measurements of spatially resolved radio jets from intermediate-mass YSOs (Anglada et al. 2018, see Table 1) and the results from high-resolution, 3D numerical simulations of MC DWs to 100 au scale (Staff et al. 2014, see Fig. 3). It is also plausible that our estimate for the jet rotational velocity is larger than the values of a few 10 km s⁻¹ typical for jets from low-mass YSOs (Bacciotti et al. 2002; Coffey et al. 2004; Ray et al. 2007). The ratio of the maser to the jet speeds depends on the density contrast between the jet and the ambient medium through V / ϒ = $\sqrt{{\rho }_{\text{jet}}/{\rho }_{\text{amb}}}$. High-angular resolution observations toward low-mass (Podio et al. 2015; Gómez-Ruiz et al. 2015) and high-mass (Anglada et al. 2018) YSOs, which are also supported by recent numerical simulations of MC-DWs (Staff et al. 2019, see Fig. 5), indicate values of $n_{{\text{H}}_2}~{10}^6$ cm⁻³ for the density of the high-velocity gas in jets. From this value and our measurement of the average ambient density in core 1 $n_{{\text{H}}_2}\approx {10}^8$ cm⁻³, we derive V / ϒ ~ 0.1. This agreeswith the ratio of the maser-averaged value of V_axi = 23 km s⁻¹ to the estimate for ϒ_axi ≈ 200 km s⁻¹.

A key prediction of the MC DW theory is that the wind mass flux is about 10% of the accretion flux (Pudritz & Ray 2019). Assuming that YSO-1 gathers most of the mass flowing toward the core cluster from the parental molecular cloud, its accretion rate can be estimated as Ṁ_acc ≈Ṁ_inf ≈ tan (i_inf) 10⁻⁴ M_⊙ yr⁻¹ (see Sect. 4.1.1). Combining the jet momentum rate Ṗ_jet ~ 3 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹ (see Sect. 4.2.1) with the jet terminal velocity ϒ_axi ≈ 200 km s⁻¹ derived above, the mass ejection rate of YSO-1 is found to be Ṁ_eje ≈ 1.5 × 10⁻⁵ M_⊙ yr⁻¹. Thus, the ratio of the ejection and the accretion rates is Ṁ_eje ∕ Ṁ_acc ≈ 0.15 cot(i_inf). For an intermediate value of (the inclination of the infall) i_inf ≈ 55°, this ratio is consistent with the value of ≈0.1 predicted by the MC DW theory.

Fig. 14 Maser positions projected perpendicular to the jet axis, Rper, vs. corresponding LOS velocities, VLOS, for the masers with VLOS ≥ −20 km s−1. Colored error bars are used to indicate values and associated errors; the colors represent the distance from YSO-1 projected along the jet axis, Raxi, as coded in the wedge on the right-hand side of the panel. The black dashed line shows the linear fit to the plotted values.

Fig. 15 Components Vaxi vs. VΦ (see definitions in Sect. 4.3) for the 3D maser velocities shown in Fig. 12. Colored error bars are used to indicate values and associated errors; the colors represent the distance from YSO-1 projected along the jet axis, Raxi, as coded in the wedge on the right-hand side of the panel. The black dashed line shows the linear fit to the subset of points within the region delimited by the horizontal and vertical black dotted lines, corresponding to VΦ ≤ 20 km s−1 and Vaxi ≥ 5 km s−1.

4.4 Change in the jet direction

Since there is no hint of multiple YSOs inside core 1 in our data, the different direction of the outflows from core 1 shown in Fig. 9 (upper panel) is most likely due to a change in the orientation of the jet ejected from YSO-1. The angle between the direction, at PA = 71°, to the far H₂ bow shock and the axis of the compact radio jet, is ≈28°. From previous studies of jet precession (see, for instance, Fendt & Zinnecker 1998; Shepherd et al. 2000) and recent 3D magnetohydrodynamics (MHD) simulations of protostar formation (see, for instance, Hirano & Machida 2019), it is known that several processes can produce a large (≥ a few 10°) change in the jet orientation: (1) radiation-induced warping of the protostellar accretion disk; (2) tidal interactions between the disk around the primary and the secondary, within a close binary system; (3) anisotropic accretion; and (4) magnetic braking. In the following, we discuss each of these points. We can readily exclude the first point, since the mass and luminosity of YSO-1 is not sufficient to warp the disk. By employing Eq. (1) of Shepherd et al. (2000), we find that the size of the observed disk is orders of magnitude smaller than the critical value beyond which the disk can be unstable to warping.

Concerning now the second point, we consider the binary system of YSOs inside the adjacent cores 1 and 3 (see Sect. 4.1.2). Assuming that the disk surface density is uniform, for Keplerian rotation, the precession frequency is given by the expression (Terquem et al. 1999, see Eq. (1)) $$w_{\text{p}}=-\frac{15}{32}\frac{M_{\text{s}}}{M_{\text{p}}}{\left(\frac{R_{\text{d}}}{D}\right)}^3\cos (\delta )\sqrt{\frac{G{M}_{\text{p}}}{R_{\text{d}}^3}},$$(4)

where M_p and M_s are the masses of the primary and secondary stars, respectively, R_d and D are the disk radius and binary separation, respectively, and δ is the inclination angle of the binary orbit with respect to the plane of the disk. From Sects. 4.1.2 and 4.2.2, we have M_p = M_⋆ ≈ 5.6 M_⊙, M_s ∕M_p ≈ 0.25, D ≈ 770 au. Referring to Fig. 10 (lower panel), R_d ≈ 250 au is estimated by taking half the spread in position (projected along the major axis of the distribution at PA = 123°) for the high-velocity molecular emissions tracing the YSO-1 disk. The line, at PA = 129°, connecting the dust emission peaks of cores 1 and 3 and the disk major axis form a small angle of ≈ 6°, which corresponds to the projection of δ on the plane of the sky: assuming cos (δ) = 1, we calculate an upper limit for w_p. Using the listed values for the parameters of Eq. (4), we infer a lower limit for the precession period T_p = 2π∕w_p ≳ 4 × 10⁵ yr. Comparing this with the dynamical timescale of the outflow t_out ~ 6.3 × 10³ yr (see Sect. 4.2.1), we could argue that is very unlikely that the jet has precessed across an angle ≥ 30° over a tinyfraction of the precession period. The precession axis should be directed close to the LOS and the binary orbit should be seen almost face-on, with the consequence that the mass of the binary system inferred from the observed V_LSR gradient would be exceedingly too high with respect to the mass of YSO-1. However, if the dynamical time underestimated the true outflow age by a factor as large as 10, the outflow lifetime could still be a significant fraction (i.e., up to one-sixth) of the precession period. Therefore, we cannot completely exclude that the change in the YSO-1 jet direction at different length scales is due to precession induced by the tidal interaction between the YSOs in cores 1 and 3.

Recent 3D nonideal MHD simulations of the collapse of a rotating, magnetized molecular clump (Tsukamoto et al. 2018; Hirano & Machida 2019; Machida et al. 2020), which aim to study the effects of a misalignment between the rotation axis of the clump and the (initially) uniform magnetic field, indicate that the evolution of the angular momentum of the central region is controlled by both anisotropic accretion and magnetic braking. While the former occurs during the earlier phases at lower density when the prestellar clump preferentially contracts along the magnetic field line due to flux freezing, the latter operates at higher density close to the protostar where an outer magnetically supported pseudo-disk and an inner centrifugally supported disk form. The above simulations show that the elongation of the disk-likestructure orbiting the protostar is approximately perpendicular to the magnetic field at large scales (10³ –10⁴ au) and to the clump rotation axis at small scales (~10² au). The direction of the protostellar outflow changes with time by a few 10° following the evolution of the accretion disk, and converges to the rotation axis of the clump once the centrifugally supported disk has grown in size.

The formation of the core cluster in IRAS 21078+5211 could have undergone processes similar to those predicted in these MHD simulations. This possibility is suggested by the finding that the most massive and evolved cores in the cluster have comparable sizes of a few 100 au and are all elongated approximately in the same direction. Indeed, Fig. 7 (upper panel) shows that core 2 is extended along an axis about parallel to the direction of both the binary system in cores 1 and 3 (at PA = 129°) and the disk around YSO-1 in core 1 (at PA = 123°). According to the aforementioned simulations, the direction, at PA ≈ 43°, of the jet from YSO-1 (the most evolved YSO in the cluster) should be close to the rotation axis of the parental clump. The large-scale magnetic field should be oriented approximately perpendicular to the SE-NW elongated, core cluster (see Fig. 1, upper panel), whose PA is 149° ± 17° (see Sect. 3.1). Therefore, in this scenario, the direction from YSO-1 to the H₂ bow shock, at PA = 71°, would be about parallel to the magnetic field. The H₂ shocks trace the direction of the jet axis at an earlier time than the radio jet because they are at a larger distance from YSO-1. The observed change in the jet orientation could then correspond to the model-predicted evolution of the protostellar outflow, from being collimated close to the magnetic field in the earliest phases to ultimately align with the rotation axis of the clump. This interpretation, relying only on marginal evidence, is speculative, but it can be tested with future measurements of the magnetic field configuration in IRAS 21078+5211.

5 Conclusions

This work combines the data of the CORE NOEMA millimeter interferometer and POETS VLBA-JVLA radio surveys to study the SF in IRAS 21078+5211 on scales from ~10⁴ to ~10 au. The CORE dust continuum and molecular line emissions allow us to map the cluster of molecular cores, identify the positions of the YSOs, and determine the kinematics of the associated accretion-ejection structures down to a few 100 au. The POETS radio continuum and water maser data complement the CORE information by unveiling the properties of the ionized gas and the 3D gas kinematics near the YSOs at radii as small as a few astronomical units. Our main results can be resumed as follows:

• The SE-NW elongated, core cluster (size ~0.1 pc) of IRAS 21078+5211 is found inside a larger (size ~1 pc) molecular cloud of similar SE-NW elongation. A V_LSR gradient of amplitude of ≈1 km s⁻¹ per 0.1 pc is detected across the major axis of the molecular cloud and the cluster. Assuming we are observing a mass flow from the harboring cloud to the cluster, we derive a mass infall rate of ≈10⁻⁴ M_⊙ yr⁻¹.

• A signature of protostellar activity, in terms of emission from high-excitation molecular lines or a molecular outflow, is found only for the most massive cores (labeled 1, 2 and 3) at the center of the cluster. The masses of the YSOs inside these three cores are estimated in the range 1–6 M_⊙. If a relevant fraction of the mass infalling onto the cluster is accreted by these more massive YSOs, it is likely that the IRAS 21078+5211 core cluster can be the progenitor of a stellar system whose most massive members are early B or late O-type high-mass stars.

• We detect a SW-NE collimated, bipolar molecular outflow emerging from the most massive core 1. It is about parallel to the double-lobe (separation ≈ 0.′′5) radio jetdetected in the POETS survey, ejected from a YSO embedded in core 1. We refer to this YSO as YSO-1. We showthat the momentum rates of the radio jet and the molecular outflow are comparable, which, also because of the similar orientation, leads us to think that the molecular outflow is driven by the radio jet.

• A V_LSR gradient of amplitude of ≈14 km s⁻¹ over 500 au,directed perpendicular to the radio jet, is revealed at the position of YSO-1. We propose an interpretation in terms of an almost edge-on, rotating disk and, by fitting a Keplerian rotation to the V_LSR pattern, we obtain a mass for YSO-1 of 5.6 ± 2.0 M_⊙.

• The water masers observed in IRAS 21078+5211 with the POETS survey are clearly tracing the YSO-1 jet. They emerge at radii of 100–300 au NE of YSO-1 and their spatial distribution and proper motions are directed within a few degrees from the axis of the radio jet.

• The radio-maser jet from YSO-1 is unique in presenting a lobe with nonthermal emission and a signature of jet rotation. The latter is the monotonic change of the maser V_LSR along a (SE-NW) direction perpendicular to the (SW-NE) collimation axis of the maser proper motions.

• We show that the maser 3D velocity pattern is consistent with the MC DW models predicting the recollimation of the wind within a rotating cylindrical flow. The masers could trace shocked gas at the wall of the cylindrical jet. We determine the jet radius to be ≈16 au and the corresponding launching radius and terminal (poloidal) velocity ≈2.2 au and ≈200 km s⁻¹, respectively. Assuming that YSO-1 gathers most of the mass flowing toward the core cluster from the parental molecular cloud, the ratio of the mass ejection and accretion rates of YSO-1 can be consistent with the value (≈0.1) predicted by the MC DW theory.

Acknowledgements

A.P. acknowledges financial support from CONACyT and UNAM-PAPIIT IN113119 grant, México. A.S.M. acknowledges support from the Collaborative Research Centre 956 (subproject A6), funded by the Deutsche Forschungsgemeinschaft (DFG) – project 184018867. D.S. acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: “Building a Habitable Earth” (SE 1962/6-1). F.B. acknowledges support from the project PRIN-INAF-MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation instruments”. H.B. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) via Sonderforschungsbereich (SFB) 881 “The Milky Way System” (sub-project B1). H.B., A.A. and S.S. acknowledge support from the European Research Council under the European Community’s Horizon 2020 framework program (2014-2020) via the ERC Consolidator grant ‘From Cloud to Star Formation (CSF)’ (project number 648505). R.G.M. acknowledges support from UNAM-PAPIIT project IN104319. R.K. acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1 and KU 2849/3-2. This work is based on observations carried out under project number L14AB with the IRAM NOEMA Interferometer and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporationpartners are: The University of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, The Leibniz Institute forAstrophysics Potsdam, and Heidelberg University; The Ohio State University, and The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota and University of Virginia.

Appendix A Velocity distributions of CH₃CN and HC₃N close to YSO-1

Figure A.1 shows that the V_LSR distributions of the CH₃CN and HC₃N emissions inproximity of YSO-1 are similar. For both molecules, the emission of the blue-shifted channels peak W-NW of YSO-1 (pinpointed by the compact 1.3 cm emission) and that of the red-shifted channels E-SE of YSO-1.

Fig. A.1 NOEMA data. The gray-scale map reproduces the 1.37 mm continuum emission, plotting values in the range 10–35 mJy beam−1. The colored dots give the (Gaussian-fitted) positions of the channel emission peaks for the CH3 CN JK = 12K–11K (K = 3–6) (upper panel) and HC3N J = 24–23 (lower panel) lines; the color denotes the channel VLSR. The yellow contours show the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels at 70, 80, and 90% of the peak emission of 0.50 mJy beam−1.

Appendix B Momentum of the molecular outflow

The momentum, P_out, of the molecular outflow can be calculated using the ¹³CO J = 2–1 and SO J_N = 6₅–5₄ emission lines, which at sufficiently high absolute LOS velocities (|V_LSR − V_sys|≥ 3 km s⁻¹; see Figs. 3 and 4) mainly trace the outflowing gas. In the limit of optically thin emission, we can use the following expressions (see, for instance, Goldsmith & Langer 1999, Eq. (10)): $$\begin{array}{l}{P}_{\text{out}}\cos (i_{\text{out}})=F{\sum }_{k=1}^2{\sum }_{n=n_{1,k}}^{n_{2,k}}{\sum }_{j=1}^J{I}_{n,j}\delta V|V_n-V_{\star }|d^2\delta \Omega ,\\ F=\frac{16\ln (2)}{B_{\text{max}}{B}_{\text{min}}}\frac{Z_{\text{mol}}(T_{\text{ex}})e^{E_{\text{u}}/(k_{\text{B}}{T}_{\text{ex}})}}{hc{g}_{\text{u}}{A}_{\text{ul}}}\frac{\mu m_{\text{H}}}{X_{\text{mol}}},\end{array}$$

where i_out is the inclination angle of the outflow with respect to the LOS, the index k considers separately the two, blue- (k = 1) and red-shifted (k = 2), outflow lobes, the sum in n extends over the velocity channels (n_1,k ≤ n ≤ n_2,k) of each lobe, and j runs over all the J pixels of the maps shown in Figs. 3 and 4. In Eq. (B.1), I_n,j and V_n are, respectively, the map intensity and V_LSR at pixel j for channel n, V_⋆ (≈ −6.4 km s⁻¹, see Sect. 4.2.2) is the V_LSR of YSO-1, δV (0.5 km s⁻¹) is the channel width, d the distance (1.63 kpc) and δΩ the solid angle of the pixel. In Eq. (B.2), B_max and B_min are the FWHM sizes of the map beam along the major and minor axes, respectively, T_ex the excitation temperature of the line, E_u, g_u and A_ul the energy and the statistical weight of the upper level and the spontaneous emission coefficient, respectively, of the transition, Z_mol(T_ex) the molecular partition function, X_mol the molecular abundance with respect to the H₂ molecule, μ = 2.8 the mean molecular weight, m_H the mass of the hydrogen atom, h the Planck constant, and c the speed oflight.

Table B.1 reports the parameters of the ¹³CO J = 2–1 and SO J_N = 6₅–5₄ transitions, recovered from databases of molecular spectroscopy. Assuming LTE conditions, the T_ex of the two lines can be approximated with the gas kinetic temperature. A lower limit is the value of ≈66 K, averaged over the full extent of the cluster, derived by Beuther et al. (2018, see Table 3) by fitting the IRAM 30 m H₂ CO data at an angular resolution of 11.′′8 with XCLASS; an upper limit is the maximum temperature of ≈200 K toward core 1 from the XCLASS fit of the CH₃CN J_K = 12_K–11_K (K = 0–6) lines (see Fig. 11). We then use for both lines T_ex = 100 K, which should be accurate to within a factor of 2. We verified that the error on the outflow momentum owing to this uncertainty on T_ex is less than a factor of 2. By fitting the emission of the two lines and the dust continuum with XCLASS, and assuming the same dust absorption coefficient and gas-to-dust mass ratio of Sect. 4.1.2, we estimated the abundance ratios (with respect to H₂) of the ¹³CO and SO molecules to be 2 × 10⁻⁶ and 10⁻⁸, respectively. Employing asingle line for each molecular species does not allow us to correct for optical depth effects, which can be significant in particular for the more abundant ¹³CO. The H₂ column density is derived from NOEMA-only data, so due to potential missing flux (by about a factor of 2 for IRAS 21078+5211, Beuther et al. 2018, see Table 3), it is likely a lower limit as well. Based on our experience with the CORE data, we judge that these estimates of the abundance ratios should be correct within a factor of 3.

The V_LSR ranges [V_1,k, V_2,k] (where V_1,k and V_2,k are the velocities of the channels n_1,k and n_2,k, respectively) for the blue- (k = 1) and red-shifted (k = 2), outflow lobes (see Table B.1) are selected in such a way to include only the channels where the ¹³CO J = 2–1 and SO J_N = 6₅–5₄ emissions are detected with a high signal-to-noise ratio (see Figs. 3 and 4; the channel at V_LSR = −28.1 km s⁻¹ of the SO J_N = 6₅–5₄ line is excluded because it is too noisy).

Summing over the two outflow lobes, we derive a momentum for the molecular outflow of ~17 and ~82 M_⊙ km s⁻¹ using the ¹³CO J = 2–1 and SO J_N = 6₅–5₄ line, respectively. In the following, we adopt the intermediate value of P_out cos(i_out) ~50 M_⊙ km s⁻¹, which, based on the above considerations, should be accurate to within an order of magnitude.

References

All Tables

All Figures

Fig. 1 Upper panel: NOEMA 1.37 mm continuum emission (color map) of IRAS 21078+5211. The white contours correspond to levels ranging from 3σ (σ = 0.5 mJy beam−1) to the peak intensity of 35 mJy beam−1, in steps of 7σ. The little black crosses indicate the positions of the seven strongest 1.37 mm peaks labeled as in Beuther et al. (2018, see Table 5), with label numbers ordered in intensity. The restoring beam of the 1.37 mm continuum map is shown in the lower left corner. Lower panel: Wide-field Infrared Survey Explorer (WISE) 4.6 μm image (color map, Wright et al. 2010) toward IRAS 21078+5211. The white contours represent the SCUBA 850 μm continuum (Di Francesco et al. 2008), showing 10 levels increasing logarithmically from 7σ (σ = 0.1 Jy beam−1) to 11.4 Jy beam−1. The magenta rectangle delimits the region plotted in the upper panel.

In the text

Fig. 2 IRAM 30 m data. Intensity-weighted velocity (color map) of the C18O J = 2–1 (upper panel), 13CO J = 2–1 (middle panel), and H2CO $J_{K_a,K_c}=$30,3–20,2 (lower panel) lines determined over a small velocity range, [−7.3, −3.5] km s−1, around Vsys = −6.1 km s−1. The plotted region corresponds to the area in which the velocity-integrated emission of the C18O J = 2–1 line is higherthan 44 Jy beam−1 km s−1, ≈10σ. In each panel, the black contours represent the SCUBA 850 μm continuum (Di Francesco et al. 2008), showing 10 levels increasing logarithmically from 1 to 11.4 Jy beam−1. The white rectangle delimits the region plotted in the upper panel of Fig. 1.

In the text

Fig. 3 Merged NOEMA–IRAM 30 m data. Each panel presents the emission of the 13CO J = 2–1 line (color map) in a different velocity channel. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted intensity (in square parentheses in milliJansky per beam) are reported. The red boxes and ticks identify the panels corresponding to the central velocities, that is |VLSR − Vsys| < 3 km s−1. The white crosses mark the positions of the seven strongest peaks of the 1.37 mm continuum emission. At the central velocities, the black dashed line gives the major axis of the 1.37 mm core cluster.

In the text

Fig. 4 Merged NOEMA–IRAM 30 m data. Each panel presents the emission of the SO JN = 65–54 line (color map) in a different velocity channel. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted intensity (in square parentheses in milliJansky per beam) are reported. The red boxes and ticks identify the panels corresponding to the central velocities, that is |VLSR − Vsys| < 3 km s−1. The white crosses mark the positions of the seven strongest peaks of the 1.37 mm continuum emission.

In the text

Fig. 5 Upper panel: NOEMA data. Velocity-integrated intensity (color map) of the CH3 CN JK = 123–113 line. The 1.37 mm continuum emission is represented with black contours, showing seven levels increasing logarithmically from 11 to 35 mJy beam−1. The positions of the three strongest peaks of the 1.37 mm continuum emission are indicated, using the same labels as in Beuther et al. (2018, see Table 5). The restoring beams of the CH3 CN JK = 123–113 line and 1.37 mm continuum maps are shown in the lower right and left corners, respectively. Lower panel: the blue and red contours reproduce the emission of the SO JN = 65–54 line integrated over the velocity ranges [−22,−13] and [−1, 5] km s−1, respectively. The plotted levels are from 0.03 to 2.5 Jy beam−1, in steps of 0.27 Jy beam−1, and from 0.6 to 2.2 Jy beam−1, in steps of 0.37 Jy beam−1, for the blue- and red-shifted SO emission, respectively. The magenta contours give the JVLA A-Array continuum at 5 cm (Moscadelli et al. 2016), showing levels at 40, 50, and 90% of the peak emission of 95 μJy beam−1: the magenta dashed line connects the two strongest (nearby, but resolved) 5 cm peaks. The two black dashed lines delimit the viewing angle of the blue-shifted emission of the SO JN = 65–54 line from the NE 5 cm peak (aligned in position with core 1, see Sect. 3.1). The black contours have the same meaning as in the upper panel.

In the text

	Fig. 6 Beam-averaged spectrum (black line) of the CH3CN JK = 12K–11K (K = 0–8) line emissiontoward the 1.37 mm continuum peak. The positions of the K = 0–8 lines are indicated. The fits (red histograms) of the K = 0–6 lines performed with XCLASS are also shown.
In the text

Fig. 7 NOEMA data. The color map reproduces the column density (upper panel) and the velocity (lower panel) of the CH3 CN emission determined by fitting with XCLASS the CH3CN JK = 12K–11K (K = 0–6) lines (including the CH3 13CN isotopologs), simultaneously. The plotted regions correspond to the areas in which the 1.37 mm continuum emission is higher than 13 and 16 mJy beam−1 for the column density and velocity plots, respectively. The 1.37 mm continuum emission is represented with black contours, showing six levels increasing logarithmically from 13 to 35 mJy beam−1: the positions of the three strongest 1.37 mm peaks are denoted using the same labels as in Beuther et al. (2018, see Table 5). The restoring beams of the 1.37 mm continuum and CH3 CN maps are shown in the lower left corner of the upper panel and lower right corner of the lower panel, respectively.

In the text

Fig. 8 NOEMA data. Upper panels: the color maps reproduce the emission of the CH3 CN JK = 125–115 (left) and CH2CO $J_{K_a,K_c}=$111,11–101,10 (right) lines. In the upper right corner of the panel, the channel VLSR (in kilometer per second) and the range of plotted values (in milliJansky per beam) are reported. The white contours reproduce the 1.37 mm continuum, with the same levels as in Fig. 7: the positions of the two nearby 1.37 mm peaks are denoted using the same labels as in Beuther et al. (2018, see Table 5), namely, 1 and 3 for the primary and secondary peak, respectively. Right panel: the blue segment connects the two peaks. In the two panels, the black contours delimit the integration areas used to produce the spectra toward cores 1 and 3, presented in the lower panel. Lower panel: spectra of the CH3 CN JK = 125–115 line toward core 1 (blue histogram) and the CH2CO $J_{K_a,K_c}=$111,11–101,10 line toward core 3 (red histogram). The flux scales for the CH3CN and CH2CO spectra are reported on the left and right, vertical axes, respectively. The blue and red dashed vertical lines indicate the approximateVLSR of the cores 1 and 3, respectively.

In the text

Fig. 9 Comparison of CORE, POETS and LBT observations. Upper panel: the gray-scale map reproduces the H2 2.12 μm emission observed with LBT toward IRAS 21078+5211. Referring to Fig. 5, lower panel, the blue and red contours have the same meaning, and the yellow contours replace the magenta contours for the JVLA 5 cm continuum. The magenta dashed, white dotted, and cyan dashed lines connect the main 5 cm peak with the secondary (nearby) 5 cm peak, the weak 5 cm spur to NE, and the bow shock traced by the LBT H2 2.12 μm emissionto SW, respectively. The white dashed box delimits the region plotted in the lower panel. Lower panel: adaptive-optics assisted LBT observations of the H2 2.12 μm emission (gray-scale map): the cyan dashed line has the same meaning as in the upper panel.

In the text

Fig. 10 NOEMA data. Upper panel: the gray-scale map reproduces the 1.37 mm continuum emission, plotting values in the range 10–35 mJy beam−1. The colored dots give the (Gaussian-fitted) positions of the channel emission peaks for the CH3 CN JK = 12K–11K (K = 3–6) and HC3N J = 24–23 lines; the colors denote the channel VLSR. The colored contours indicate the half-peak level for the CH3CN JK = 123-113 emission in individual channels. The yellow contours show the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels at 70, 80, and 90% of the peak emission of 0.50 mJy beam−1. Lower panel: the gray-scale map, yellow contours, and colored dots have the same meaning as in the upper panel. Only the channel emission peaks at the extreme blue- and red-shifted velocities (reported in the lower right corner of the panel) are shown. The white continuous and dashed lines indicate the best-fit major axis of the distribution of the peaks and the corresponding fit uncertainty, respectively.

In the text

Fig. 11 NOEMA data. The color map shows the distribution of the rotational temperature derived by fitting with XCLASS the emission of the CH3CN JK = 12K–11K (K = 0–6) lines (including the CH3 13CN isotopologs). The plotted region corresponds to the area in which the 1.37 mm continuum emission is higher than 16 mJy beam−1. The black contours reproduce the 1.37 mm continuum, showing the same levels as in Fig. 7: the positions of the three strongest 1.37 mm peaks are indicated using the same labels as in Beuther et al. (2018, see Table 5). The restoring beam of the CH3 CN maps is shown in the lower right corner of the panel.

In the text

Fig. 12 POETS water maser observations. Upper panel: colored dots and arrows give absolute positions and proper motions of the 22 GHz water masers determined with VLBI observations (Moscadelli et al. 2016); colors denote the maser VLSR. The dot area scales logarithmically with the maser intensity. The black contours indicate the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels from 10 to 90% of the peak emission of 0.50 mJy beam−1 in steps of 10%. Lower panel: colored dots and black contours have the same meaning as in the upper panel. Colored cones are employed to visualize the maser 3D velocities, representing the inclination with respect to the LOS through the ellipticity of the cone basis and the uncertainty in the direction by means of the cone aperture.

In the text

Fig. 13 χ2 values from the Keplerian fit (see Eq. (1)) vs. the free parameter M⋆ sin 2(irot). The horizontal black dashed line indicates the value ${\chi }_{\text{min}}^2$+1.0 = 4.12, which corresponds to the 1σ confidence level for the χ2-distribution with onefree parameter (Lampton et al. 1976). The two vertical blue dashed lines indicate the 1σ confidence interval for the best-fit value of the free parameter.

In the text

Fig. 14 Maser positions projected perpendicular to the jet axis, Rper, vs. corresponding LOS velocities, VLOS, for the masers with VLOS ≥ −20 km s−1. Colored error bars are used to indicate values and associated errors; the colors represent the distance from YSO-1 projected along the jet axis, Raxi, as coded in the wedge on the right-hand side of the panel. The black dashed line shows the linear fit to the plotted values.

In the text

Fig. 15 Components Vaxi vs. VΦ (see definitions in Sect. 4.3) for the 3D maser velocities shown in Fig. 12. Colored error bars are used to indicate values and associated errors; the colors represent the distance from YSO-1 projected along the jet axis, Raxi, as coded in the wedge on the right-hand side of the panel. The black dashed line shows the linear fit to the subset of points within the region delimited by the horizontal and vertical black dotted lines, corresponding to VΦ ≤ 20 km s−1 and Vaxi ≥ 5 km s−1.

In the text

Fig. A.1 NOEMA data. The gray-scale map reproduces the 1.37 mm continuum emission, plotting values in the range 10–35 mJy beam−1. The colored dots give the (Gaussian-fitted) positions of the channel emission peaks for the CH3 CN JK = 12K–11K (K = 3–6) (upper panel) and HC3N J = 24–23 (lower panel) lines; the color denotes the channel VLSR. The yellow contours show the JVLA A-Array continuum at 1.3 cm (Moscadelli et al. 2016), showing levels at 70, 80, and 90% of the peak emission of 0.50 mJy beam−1.

In the text