IRAS23385+6053: An embedded massive cluster in the making

TL;DR

This study investigates IRAS23385+6053, a high-mass star-forming region, revealing multiple cores, a potential circumstellar disk, and signs of ongoing collapse, indicating an early stage of massive cluster formation.

Contribution

First high-resolution multi-line observations of IRAS23385+6053, identifying multiple cores, a possible circumstellar disk, and dynamic motions in a very early massive star formation stage.

Findings

Six distinct cores identified in the clump.

Detection of a potential circumstellar self-gravitating disk.

Evidence of a bipolar outflow driven by a low-mass star.

Abstract

This study is part of the project ``CORE'', an IRAM/NOEMA large program consisting of observations of the millimeter continuum and molecular line emission towards 20 selected high-mass star forming regions. We focus on IRAS23385+6053, which is believed to be the least evolved source of the CORE sample. The observations were performed at ~1.4 mm and employed three configurations of NOEMA and additional single-dish maps, merged with the interferometric data to recover the extended emission. Our correlator setup covered a number of lines from well-known hot core tracers and a few outflow tracers. The angular (~0.45"$-$0.9") and spectral (0.5 km/s) resolutions were sufficient to resolve the clump in IRAS23385+6053 and investigate the existence of large-scale motions due to rotation, infall, or expansion. We find that the clump splits into six distinct cores when observed at sub-arcsecond resolution. These are identified through their 1.4 mm continuum and molecular line emission. We produce maps of the velocity, line width, and rotational temperature from the methanol and methyl cyanide lines, which allow us to investigate the cores and reveal a velocity and temperature gradient in the most massive core. We also find evidence of a bipolar outflow, possibly powered by a low-mass star. We present the tentative detection of a circumstellar self-gravitating disk lying in the most massive core and powering a large-scale outflow previously known in the literature. In our scenario, the star powering the flow is responsible for most of the luminosity of IRAS23385+6053 (~$3000~L_\odot$). The other cores, albeit with masses below the corresponding virial masses, appear to be accreting material from their molecular surroundings and are possibly collapsing or on the verge of collapse. We conclude that we are observing a sample of star-forming cores that is bound to turn into a cluster of massive stars.

Tables

Table 1. Table 1: List of tracers detected in each core and corresponding parameters. The question marks indicate uncertain (non)detections, mostly due to the lack of a clear peak of emission coinciding with the core.

tracer	$\nu$ (MHz)	$E_{\mathrm{up}}$ (K)	$S{\mu }^2$ (D2)	A1	A2	B	C	D	E
H2CO($3_{03}$–$2_{02}$)	218222.192	20.96	16.3	Y	Y	Y	Y	Y	Y?
H2CO($3_{22}$–$2_{21}$)	218475.632	68.09	9.06	Y	Y	Y	Y	Y	Y?
CH3OH($4_2$–$3_1$) E1	218440.050	34.98	3.48	Y	Y	Y	Y	Y	Y?
HC3N(24–23)	218324.723	130.98	334	Y	Y?	Y	N	Y	Y?
OCS(18–17)	218903.356	99.81	9.21	Y	N?	N?	N	N	N
CH3CN(12–11)	220747.261a	68.87a	254a	Y	Y?	N?	N	N	N?
H2CCO(${11}_{1,11}$–${10}_{1,10}$)	220177.569	76.46	66.0	Y	N?	N	Y	N	N
1.4 mm continuum	218917	—	—	Y	Y	Y	N?	Y?	Y

Table 2. Table 2: Parameters of the cores in the IRAS 23385+6053  region.

core	$\Delta \alpha$a	$\Delta \delta$a	$R_c$	$V_{\mathrm{LSR}}$b	$\Delta V$b	$T_{\mathrm{rot}}$	$S_{\nu }$	$M_{\mathrm{vir}}$c	$M_{H_2}$	$n_{H_2}$	$\frac{M_{H_2}}{M_{\mathrm{vir}}}$	${\dot{M}}_{\mathrm{acc}}$	$t_{\mathrm{acc}}$
	(″)	(″)	(″)	(km s-1)	(km s-1)	(K)	(mJy)	($M_{\odot }$)	($M_{\odot }$)	(cm-3)		($M_{\odot }$ yr-1)	(yr)
A1	0.81	0.00	0.69	–50.6	3.2	39	62	35	24	$1.8\times {10}^7$	0.67	$3.1\times {10}^{-3}$	$7.6\times {10}^3$
A2	0.28	–0.34	0.43	–50.4	2.8	33	25	17	12	$3.8\times {10}^7$	0.70	$2.1\times {10}^{-3}$	$5.7\times {10}^3$
B	2.36	0.44	0.81	–51.6	1.9	31	36	15	18	$9.0\times {10}^6$	1.25	$6.5\times {10}^{-4}$	$2.8\times {10}^4$
C	–0.46	–0.79	0.43	–50.8	2.7	27	3.8	16	2.2	$7.3\times {10}^6$	0.14	$1.9\times {10}^{-3}$	$1.2\times {10}^3$
D	0.69	1.38	0.69	–51.3	3.6	26	12	45	7.4	$5.7\times {10}^6$	0.17	$4.4\times {10}^{-3}$	$1.7\times {10}^3$
E	–0.59	0.15	0.47	–50.2	2.0	24	10	9.4	6.8	$1.7\times {10}^7$	0.72	$7.6\times {10}^{-4}$	$9.0\times {10}^3$

Equations

$$V(x,y)=V_{\mathrm{sys}}+{\text{varv}}_{\mathrm{rot}}{\left(\frac{R}{R_{\mathrm{o}}}\right)}^{\alpha }\frac{x}{R}\sin \theta +{\text{varv}}_{\mathrm{rad}}{\left(\frac{R}{R_{\mathrm{o}}}\right)}^{\beta }\frac{y}{R}\tan \theta$$

$$P\cos i=\sum _j\frac{N_j}{X_{\mathrm{C}{\mathrm{H}}_3\mathrm{CN}}}(V_j-V_{\mathrm{sys}})\mu m_{\mathrm{H}}{d}^2\delta \Omega$$

$$M_{\mathrm{vir}}=\frac{5}{6}\frac{R_{\mathrm{c}}}{G}\frac{(\Delta V{)}^2}{8\ln 2}\left(1+\sqrt{1-\frac{p_{\mathrm{e}}}{p_0}}\right)$$

Full text

*Aims. * We focus on IRAS 23385+6053, which is believed to be the least evolved source of the CORE sample. This object is characterized by a compact molecular clump that is IR dark shortward of 24 $\mu$m and is surrounded by a stellar cluster detected in the near-IR. Our aim is to study the structure and velocity field of the clump.

Methods. * The observations were performed at $\sim$1.4 mm and employed three configurations of NOEMA and additional single-dish maps, merged with the interferometric data to recover the extended emission. Our correlator setup covered a number of lines from well-known hot core tracers and a few outflow tracers. The angular ($\sim$0$\aas@@fstack{\prime\prime}$45–0$\aas@@fstack{\prime\prime}$9) and spectral (0.5 km s-1*) resolutions were sufficient to resolve the clump in IRAS 23385+6053 and investigate the existence of large-scale motions due to rotation, infall, or expansion.

*Results. * We find that the clump splits into six distinct cores when observed at sub-arcsecond resolution. These are identified through their 1.4 mm continuum and molecular line emission. We produce maps of the velocity, line width, and rotational temperature from the methanol and methyl cyanide lines, which allow us to investigate the cores and reveal a velocity and temperature gradient in the most massive core. We also find evidence of a bipolar outflow, possibly powered by a low-mass star.

*Conclusions. * We present the tentative detection of a circumstellar self-gravitating disk lying in the most massive core and powering a large-scale outflow previously known in the literature. In our scenario, the star powering the flow is responsible for most of the luminosity of IRAS 23385+6053 ($\sim$3000 $L_{\odot}$). The other cores, albeit with masses below the corresponding virial masses, appear to be accreting material from their molecular surroundings and are possibly collapsing or on the verge of collapse. We conclude that we are observing a sample of star-forming cores that is bound to turn into a cluster of massive stars.

††offprints: R. Cesaroni,

1 Introduction

Studies of high-mass ($\ga 10^{4}\leavevmode\nobreak\ L_{\odot}$) star formation have been hindered until recently by limited angular resolution and sensitivity. Young early-type stars are born deeply embedded in their parental cocoons, which makes observations longward of a few 100 $\mu$m necessary to investigate the gas and dust emission from their natal environment. In turn, this makes interferometric observations at (sub)mm wavelengths the only possibility to achieve the sub-arcsecond resolutions needed to study objects typically located at a distance of several kpc.

The advent of instruments such as the Atacama Large Millimeter/submillimeter Array (ALMA), the NOrthern Extended Millimeter Array (NOEMA), and the upgraded Submillimeter Array (SMA) has allowed a breakthrough in this respect. Not only do the superior angular resolution and sensitivity make it possible to resolve structures $\la$200 au up to distances of several kpc, but the large number of antennas allows good sampling of the $u,\varv$ plane in a short time, enabling surveys of a large number of targets.

With this in mind, we undertook the large program “CORE” (P.I. Henrik Beuther) with the IRAM interferometer, NOEMA, targeting 20 high-mass star-forming regions. There are multiple goals of the project: to search for circumstellar disks and study their properties; to establish the level of fragmentation of the pc-scale molecular clumps harbouring the massive young stellar objects (YSOs); and to understand the chemical composition of the gas. While some general results for the whole sample have already been presented in a previous paper (Beuther et al. beu18 (2018); hereafter BEU18), in the present study we focus on the particular source IRAS 23385+6053, in the footsteps of other articles of the same sample (Ahmadi et al. ahmadi2018 (2018); Mottram et al. submitted; Bosco et al. submitted).

IRAS 23385+6053, located at a kinematic distance of 4.9 kpc (Molinari et al. mol98 (1998)), is also known as Mol160 and was selected as a candidate massive protostar by Palla et al. (palla91 (1991)), on the basis of its IRAS colours. Association with water maser emission (Palla et al. palla91 (1991)) revealed star formation activity, later confirmed by other studies such as Molinari et al. (mol96 (1996)), who detected ammonia emission from the source, indicating the presence of dense molecular gas associated with it. Subsequently, Molinari et al. (mol98 (1998)) investigated the structure of the source by comparing interferometric maps of the continuum and molecular line emission with images of the continuum emission in the mid-IR. From this comparison one sees that the high-mass object is embedded in an IR-dark core surrounded by an IR-luminous stellar cluster (see Fig. A.26 of Faustini et al. fau09 (2009)). While the IRAS luminosity of the region is \sim$$1.6\times 10^{4}\leavevmode\nobreak\ L_{\odot}, MIPSGAL/Spitzer data (Molinari et al. 2008b ) have demonstrated that the contribution of the core amounts to only \sim$$3\times 10^{3} $L_{\odot}$, consistent with IRAS 23385+6053 being a massive protostar in the main accretion phase (Molinari et al 2008a ). The rest of the luminosity arises from the surrounding cluster, which is also likely ionizing the two Hii regions detected by Molinari et al. (mol02 (2002)). Figure 1a shows a pc-scale view of the region with the molecular core, the Hii regions, and the far-IR emission tracing the stellar cluster.

The presence of a compact molecular core was confirmed by the interferometric observations at 3 and 1.3 mm of Fontani et al. (fonta04 (2004)) and, more recently, Wolf-Chase et al. (wolf12 (2012)). From their observations it was also possible to identify a main core and a secondary, less prominent core separated by $\sim$2″ to the NE. The temperature of the main core was estimated to be quite low ($\sim$40 K) compared to the typical values of hot molecular cores (100–200 K; see Kurtz et al. ppiv (2000), Cesaroni cesa05b (2005)) and it was concluded that IRAS 23385+6053 is indeed a rapidly accreting massive YSO in a very early stage of its evolution. Since accretion is tightly connected to ejection, this conclusion appears confirmed by the detection of outflow tracers such as broad line wings and class I methanol maser emission associated with the cores (Molinari et al. mol98 (1998); Zhang et al. zha01 (2001, 2005); Kurtz et al. kur04 (2004); Wolf-Chase et al. wolf12 (2012)), although a clear bipolar outflow pattern (oriented NE–SW) is seen only on large scales ($\sim$1$\aas@@fstack{\prime}$5) in the single-dish maps of Wu et al. (wu05 (2005)). This orientation is also roughly consistent with the class I CH3OH masers, which are distributed over a region of 6″$\times$3″elongated in the NNE–SSW direction and centered on the molecular core detected by Molinari et al. (mol98 (1998)). In summary, all these characteristics make IRAS 23385+6053 in all likelihood the youngest object of the CORE sample.

In the present study, we make use of the new high-resolution and high-sensitivity NOEMA observations to pursue the investigation of IRAS 23385+6053 and establish its structure on scales of $\sim$0.02 pc. After describing the observations in Sect. 2, we present the results in Sect. 3, while Sect. 4 is devoted to the derivation of the physical parameters of the cores and a detailed analysis of each of them. In Sect. 5 we discuss a possible scenario for IRAS 23385+6053 and compare its properties to those of similar regions. Finally, a summary is given in Sect. 6.

2 Observations

2.1 IRAM/NOEMA interferometer

The region IRAS 23385+6053 was observed in the framework of the NOEMA large program CORE. Details about the program and sample can be found in BEU18. Sources were always observed in pairs of two in the track-sharing mode. In this case, IRAS 23385+6053 was observed together with G108.7575–0.986 in three configurations (A, C and D) between January 2015 and October 2016.  The number of antennas in the array varied between 6 and 8. The phase center for IRAS 23385+6053 was $\alpha$(J2000.0) 23h40m54$\aas@@fstack{s}$400  and $\delta$(J2000.0) +61°10′28$\aas@@fstack{\prime\prime}$020 with the velocity of rest of –50.2 km s*-1*. As gain calibrator, we used the quasar 0059+581, and bandpass and flux calibrators were 3C 454.3 and MWC 349.

The spectral range covered with the broad-band WIDEX correlator was between 217.167 and 220.834 GHz at 1.95 MHz spectral resolution (corresponding to $\sim$2.7 km s*-1*). Furthermore, eight narrow-band high-spectral resolution units ($\sim$0.5 km s*-1*) were distributed over the bandpass, mainly covering lines from CH3CN and H2CO (for more details see Table 2 in BEU18).

Data calibration and imaging were mainly conducted within the GILDAS framework, with the calibration package CLIC, and the imaging program MAPPING. The continuum was produced by collapsing the line-free part of the spectrum into a single continuum channel. Furthermore, the continuum data were self-calibrated in CASA to improve the signal-to-noise ratio. The continuum $1\sigma$ rms of the final images is 0.11 mJy beam*-1*, and for the line cubes we achieved a typical rms of $\sim$6 mJy beam*-1* in a 0.5 km s*-1* channel. The spatial resolutions of the continuum and line data are 0$\aas@@fstack{\prime\prime}$46 and $\sim$0$\aas@@fstack{\prime\prime}$9, respectively, due to the different weighting adopted (uniform for the continuum, robust=5 for the lines). The absolute flux scale is estimated to be correct to within 20%.

2.2 IRAM/30-m telescope

The above interferometer data were complemented with short spacing observations at the IRAM 30-m telescope. Details about the 30-m observations can be found in Ahmadi et al. (ahmadi2018 (2018)) and Mottram et al. (submitted). Each source was observed in the 1.4 mm band in the On-The-Fly (OTF) observation mode in map sizes of typically 1′$\times$1′. OTF maps were conducted in both right ascension and declination directions to reduce scanning effects. Calibration of the single-dish data was done in CLASS, and the merging of the single-dish and interferometer data was again conducted within the MAPPING program. While most compact structures were still cleaned with the CLARK-algorithm, the more extended 13CO(2–1) data were imaged with the Steer-Dewdney-Ito (SDI) method. Details about the merging and imaging are given in Mottram et al. (submitted). The resulting angular resolution and $1\sigma$ rms in a 3 km s*-1* channel are $\sim$0$\aas@@fstack{\prime\prime}$46 and 2 mJy beam*-1*, respectively.

3 Results

Below we present our main findings. Unless otherwise specified, the maps used in our study are those obtained with the NOEMA interferometer, because only the emission of the CO isotopologues may require merging the NOEMA and 30-m data to be properly imaged.

3.1 Continuum emission

A study of the continuum emission from IRAS 23385+6053 and the other sources of the CORE sample has already been presented by BEU18. Using the clumpfind algorithm, they identified 3 cores, whose total flux density is $\sim$180 mJy. This is greater than the value of $\sim$150 mJy obtained by Fontani et al. (fonta04 (2004)). Such a difference is larger than the noise of Fontani et al.’s map ($\sim$2 mJy/beam) and can be due to the different $u,\varv$ coverage of the two data sets. In fact at the time of Fontani et al.’s observations the Plateau de Bure interferometer was equipped with only 5 antennas and 2 configurations of the array were used, compared with our NOEMA observations performed with a number of antennas ranging from 6 to 8 and 3 different configurations. In particular, unlike Fontani et al., we used also the most compact (D-array) configuration and it is thus not surprising that part of the flux density present in our images is filtered out in their data. We show in Fig. 2 the map of the continuum emission, with overlaid the positions of the three cores identified by BEU18.

3.2 Line emission

While the correlator setup covers many potentially useful molecular lines, only some of these have been detected with a typical 1$\sigma$ noise of 3 mJy/beam ($\sim$0.1 K) in a 3 km s*-1* channel. The most prominent are the 13CO and C18O (2–1) transitions. However, these CO isotopologues are largely sensitive to extended emission, which masks the compact emission from the core. In Fig. 1 we compare a map of the mean emission in the 13CO(2–1) line to the Spitzer/MIPS image of the region at 24 $\mu$m. More precisely, in Fig. 1a we show the 13CO map made with the 30-m telescope, while in Fig. 1b the same map has been obtained after merging the 30-m and NOEMA data. One can see that the compact core is embedded in a larger emitting region extending over at least $\sim$20″ (0.48 pc). As noted by Molinari et al (2008b ), the core is detected at 24 $\mu$m (whereas it is IR-dark at shorter wavelengths). It is worth noting that no clear bipolar outflow structure can be identified from our CO data, neither in the channel maps nor in the spectra. Although this might be due to the limited spectral resolution (3 km s*-1*) and the relatively small region covered by the observations, one must also consider that 12CO would be a better outflow tracer than its isotopologues and the outflow orientation could also prevent detection of it, if the axis lies close to the plane of the sky and/or the gas is ejected at low velocity. The 13CO emitting region appears to lie between the two compact Hii regions (dotted contours in Fig. 1) imaged at 3.6 cm by Molinari et al. (mol02 (2002)), which are ionized by B-type stars belonging to the surrounding cluster.

The present study is focused on compact cores. These are much better traced by other molecular species than the CO isotopologues. Our correlator setup covers several transitions of different molecular species, as one can see in the broad-band spectrum shown in Fig. 3, obtained with the WIDEX correlator towards the peak of the emission. However, for our purposes we will focus mostly on the lines covered also by the high spectral resolution (0.5 km s*-1*) units of the narrow-band correlator and marked in red in Fig. 3. These belong to formaldehyde (H2CO), methanol (CH3OH), methyl cyanide (CH3CN), cyanoacetylene (HC3N), carbonyl sulfide (OCS), and ketene (H2CCO). The maps of the emission averaged over the lines of these species are shown in Fig. 4. Visual inspection of these images and the continuum emission map in Fig. 2, allows us to identify 6 cores, whose positions are marked in Fig. 4 with crosses.

More in detail, cores A1+A2, B, and E have been found by BEU18 from the continuum map and are clearly seen in the CH3OH map in Fig. 4, core D is well defined in the CH3OH map, and core C appears as a separate entity only in the H2CCO map. We point out that we have decided to split core 1 in Table A.1 of BEU18 into two cores, A1 and A2, because the head-tail structure seen in the continuum is not reproduced in some of the lines. In particular, the HC3N, OCS, and CH3CN emission seems to trace a circular, barely resolved core located at the position of the continuum peak (A1).

Not all molecules are detected in all cores. Table 1 lists the tracers revealed in each core, while Fig. 5 illustrates the locations of the cores through a sketch, to facilitate core identification and comparison with the maps.

An interesting feature that can be seen in the maps of the H2CO and CH3OH emission is the presence of two blobs, F1 and F2 (see Fig. 6), separated by $\sim$10″ or 0.24 pc and symmetrically disposed to the SE and NW of core C. This is shown in Fig. 6, where the dashed line joins the two blobs. While these might be two additional cores, the non detection of continuum emission from them hints at a different nature. We propose that they could be the tips of the lobes of a bipolar outflow originating from C or a nearby, low-mass undetected core. Indeed, the sensitivity (at a 5$\sigma$ level) of our continuum map is $\sim$0.1–0.5 $M_{\odot}$ for a dust temperature in the range 20–70 K. These values are typically above the mass of envelopes/disks around class I YSOs. We will further discuss this outflow hypothesis in Sect. 4, whereas in Sect. 4.2 we will comment on the possible existence of another outflow powered by a YSO in core A1.

4 Analysis

The goal of our study is to establish the nature of the cores identified in our maps. In the following we derive the core physical parameters and analyse each of them in more detail.

4.1 Physical parameters of the cores

The observed ratio between the brightness temperatures of the H2CO ($3_{22}$–$2_{21}$) and ($3_{03}$–$2_{02}$) transitions over the cores is in all cases $>$0.56, the maximum LTE ratio expected in the optically thin case111 If the emission is thin, the ratio between the brightness temperatures of two lines is equal to the ratio between the corresponding optical depths, therefore: $T_{\rm B}^{2}/T_{\rm B}^{1}\simeq[(S\mu^{2})_{2}/(S\mu^{2})_{1})]\,(\nu_{2}/\nu_{1})\,\exp[-(E_{2}-E_{1})/T]<[(S\mu^{2})_{2}/(S\mu^{2})_{1})]\,(\nu_{2}/\nu_{1})\simeq 0.56$, where the line strengths and frequencies are given in Table 1. . This indicates that the lines must be optically thick and their ratio cannot be used to estimate the gas temperature and column density of the cores. For the same reason, H2CO is not suitable for the study of the velocity field, because it traces the envelope around the cores as suggested by the size of the H2CO emission in Fig. 4. From this figure, one sees that most of the cores are recognizable in the methanol map and we thus prefer to use this species to obtain a picture of the core physical properties. In Fig. 8 we show maps of the line velocity and full width at half maximum obtained from the 1st and 2nd moments of the CH3OH line. Overlaid on this is the map of the zero moment (integrated intensity) in the same line.

As one can see from Fig. 3, two more lines of methanol are detected in the bandwidth covered by the WIDEX correlator: the $8_{0}$–$7_{1}$ E2 and $20_{1}$–$20_{0}$ E1 transitions. One can use these transitions to calculate the rotational temperature, $T_{\rm rot}$, of the CH3OH gas all over the cores. For this purpose we have used the eXtended CASA Line Analysis Software Suite (XCLASS) tool222 XCLASS is available at https://xclass.astro.uni-koeln.de

(Möller et al. xclass (2017)), which simultaneously fits the lines of a molecular species assuming LTE by varying the relevant physical parameters, namely the source angular size, rotational temperature, column density of the molecule, systemic LSR velocity, and line width. One of the advantages with respect to rotation diagrams is that XCLASS takes into account the line opacities in the fit. In our case, the region over which CH3OH is detected is much greater than the synthesized beam and we have thus assumed a beam filling factor equal to 1. The maps of $T_{\rm rot}$ and CH3OH column density are shown in Fig. 7, where both quantities peak towards the position of A1.

The same method can be applied to the CH3CN(12–11) $K$=0–6 transitions, which are known to be excellent hot core tracers. As such, these are well suited to estimate the physical parameters of chemically rich cores like A1 and A2 and indeed the CH3CN emission is detected only towards these two cores. We present the corresponding $T_{\rm rot}$ and column density maps in Fig. 9. The obvious difference with respect to Fig. 7 is that the temperature obtained from CH3CN is about twice as much as that derived from CH3OH. While this could be due to methanol being sub-thermally excited (Bachiller et al. bachi95 (1995); Kalenskii & Kurtz kaku16 (2016)), it is most likely that the interplay between opacity and temperature gradients plays a dominant role. Unlike the other cores, A1 coincides with the peak of column density and it is thus possible that species with different abundances (and opacities) trace regions with different temperatures.

Now we compute the mean values of the physical parameters of the cores. For this purpose, it is necessary to establish the border of each core. In order to simplify the problem, we assume the cores to be spherical, so that all we need to estimate is the core radius. In Table 2 we give the peak positions of the cores determined from the line and continuum maps. In particular, the positions of cores A1, B, and E are taken from Table A.1 of BEU18. These authors used clumpfind to determine the core sizes, but this method was applied to the continuum image where only 3 cores were identified. Here we prefer to use a different approach. For each core, we consider the separation from the nearest of the other cores and assume that this is twice the radius of it. We give the radii, $R_{\rm c}$, in Table 2. These values should be considered with two caveats in mind. On the one hand, the method used by us overestimates the radius, because two neighbourhing cores are not necessarily in contact. On the other hand, we underestimate the separation because we see only the projection of it on the plane of the sky.

We calculate the mean values of the LSR velocity ($V_{\rm LSR}$), line full width at half maximum ($\Delta V$), and temperature ($T_{\rm rot}$) by averaging the corresponding CH3OH parameters (see Figs. 7 and 8) over the surface of each core. Finally, we obtain the total core flux densities ($S_{\nu}$) by integrating the continuum emission (Fig. 2) inside the borders of the cores. All these quantities are listed in Table 2, where we give also virial masses, $M_{\rm vir}$, and the core masses, $M_{\rm H_{2}}$, computed from the temperature, line widths, and flux densities according to Eq. (3) of MacLaren et al. (mclar (1988)) and Eq. (1) of Schuller et al. schul09 (2009), under the assumptions of constant density, gas-to-dust mass ratio of 150, and dust opacity of 0.9 cm2 g*-1* at 1.4 mm. In the same table we report the H2 densities, for a mean molecular weight of 2.8, and the ratios $M_{\rm H_{2}}/M_{\rm vir}$.

Clearly, the velocity changes significantly from core to core, with cores A1, A2, and E being red-shifted by at least $\sim$1 km s*-1* with respect to B and D, while cores A1 and D are those with the largest line width and possibly the highest level of turbulence. We note that in all cases the line is wider than expected from pure thermal broadening ($<$0.42 km s*-1* for $T<120$ K – see Fig. 9), which implies a contribution from non-thermal motions. The line width is also large close to core C (see Fig. 8b), where a blue-shifted “spot” is seen in the velocity map (Fig. 8a; see also Sect. 4.5). Whether core C is really associated with such a feature is questionable. We remind the reader that this core has been identified only from the faint H2CCO emission, which makes it difficult to determine a precise position for it. In any case, it appears that close to the geometrical center of the (putative) outflow there is a small region with a large velocity dispersion, which could contain the source powering the flow. An alternative possibility is that this “spot” is where the flow from C impinges on the dense gas, thus causing the observed enhancement of the line width. It is worth noting that Wolf-Chase et al. (wolf12 (2012)) have detected an H2 knot (MHO 2921 in their notation) close to core C, which is probably due to emission from shocks.

The velocity difference between blobs F1 and F2 is small, only $\sim$1 km s*-1*, but this is expected if the outflow lobes lie close to the plane of the sky. However, very little velocity dispersion (the line FWHM is $\sim$2 km s*-1*) is observed towards these blobs, in contrast with the expected line broadening at the tips of a bipolar flow. We will further discuss the outflow hypothesis in Sect. 4.5.

4.2 Core A1

As already explained, we prefer to identify two cores (A1 and A2) where the analysis of BEU18 finds only one. This choice, based on the morphology of the line emission, is further confirmed by the different kinematical properties, recognizable in Fig. 8 especially for the line FWHM which is larger (by $\sim$1 km s*-1*) in A1, compared to A2.

The velocity field in A1 can be better analysed by means of high density/temperature tracers such as CH3CN. In Fig. 10 we show a map of the peak velocity of the CH3CN(12–11) transition obtained by fitting Gaussian profiles simultaneously to the strongest ($K=0$ to 4) components, after fixing the line separations to the laboratory values and forcing the line widths to be identical. This method has the advantage of reducing the velocity and line width uncertainties with respect to fitting each $K$ component independently or computing the first and second moments of the lines. The velocity pattern presents two red-shifted and two blue-shifted peaks defining two velocity gradients along directions (see the dotted and dashed lines in Fig. 10) roughly perpendicular to each other. The CH3CN emission (see contour map in Fig. 10) clearly peaks towards A1, which proves that this emission is largely dominated by A1 and the CH3CN velocity is affected only marginally by the gas in A2. In the following, we discuss three scenarios that could explain the observed velocity field in core A1.

4.2.1 Keplerian disk with expansion

The presence of two velocity gradients, directed SE–NW and SW–NE, could be due to the combination of Keplerian rotation and expansion along the surface of the disk, a situation reminiscent of disk winds.

It is worth noting that we cannot distinguish between expansion and infall on the basis of the observed velocity pattern, since in both cases the projected velocities along the line of sight give origin to the same type of velocity field. In practice, reversing the sign of the radial velocity from positive (for expansion) to negative (for infall) only swaps the red- and blue-shifted sides of the disk, which can be compensated with a rotation by 180° of the position angle of the disk. Here, we rule out infall a priori, because the existence of a pair of blue- and red-shifted velocity peaks close to the disk border can be explained only if the radial velocity component increases with radius. This situation is in contrast with infall, which accelerates towards the star, whereas it is acceptable for expansion, where the material is accelerated outward.

While a physical model to reproduce the observed velocity and intensity goes beyond the scope of the present article, we can compare the velocity map in Fig. 10 with the map of the line-of-sight velocity computed from a purely kinematical model assuming Keplerian rotation and radial expansion in a geometrically thin disk. The observed velocity is given by the expression

[TABLE]

where $x$ and $y$ are cartesian coordinates in the plane of the sky with $y$ along the the projected major axis of the disk, $V_{\rm sys}$ is the systemic velocity with respect to the local standard of rest (LSR), $R=\sqrt{x^{2}+(y/\cos\theta)^{2}}$, $R_{\rm o}$ is the outer radius of the disk, $\varv_{\rm rot}$ and $\varv_{\rm rad}$ are the azimuthal and radial velocity components in the plane of the disk at $R_{\rm o}$, $\theta$ is the angle between the disk axis and the line of sight, and we have assumed $\alpha=-1/2$ for Keplerian rotation and $\beta=1$ for expansion. In Fig. 11 we show the velocity map obtained for fiducial values of the parameters ($\theta=30\degr$, $\varv_{\rm rot}=0.7$ km s*-1*, $\varv_{\rm rad}=2$ km s*-1*, $V_{\rm sys}=-50.7$ km s*-1*), obtained by visual comparison with Fig. 10. In particular, the inclination $\theta$ has been obtained from the ratio between the minor axis ($\sim$2″) and the major axis ($\sim$2$\aas@@fstack{\prime\prime}$3) of the velocity map in Fig. 10, as $\theta=\arccos(2/2.3)$. Note that all lengths are normalized with respect to $R_{\rm o}$ and the velocity is computed down to a minimum radius $R_{\rm i}=0.15\,R_{\rm o}$.

The model presents a number of features which are clearly not in agreement with the data. For example, the two inner velocity peaks are too close to the center of the disk and the directions of the two velocity gradients (dashed and dotted lines in Fig. 11) are not perpendicular to each other as they are in Fig. 10. However, we stress that the comparison between Fig. 10 and Fig. 11 is only qualitative, as our calculations do not take into account the angular resolution of the observations and no attempt is made to compute the line intensity and width. We do not want to fit the data but only to show that rotation plus expansion could mimic the existence of two pairs of blue- and red-shifted velocity peaks (connected by the dotted and dashed lines in Fig. 11), one located inside the disk, the other close to the border of it. These double peaks generate the twist in the velocity pattern that is seen in the data.

Since in this model the disk is undergoing Keplerian rotation, one can compute the stellar mass from the high-velocity peaks close to the center of the disk. From Fig. 10 one obtains a rotation velocity333 The velocity is obtained as half the difference between the velocities observed at the blue-shifted ($\sim$–51 km s*-1*) and red-shifted ($\sim$–50 km s*-1*) peaks located along the dotted line in Fig. 10.

of 0.5 km s*-1* at a radius of 0$\aas@@fstack{\prime\prime}$73 (or 3600 au), which implies a dynamical mass \mbox{M_{\rm dyn}}\sim1 $M_{\odot}$. This value appears by far too small compared to the mass of the core ($\sim$24 $M_{\odot}$; see Table 2), to satisfy the condition $M_{\rm star}>M_{\rm disk}$ for Keplerian rotation. Moreover, it is very unlikely that a solar type star associated with a chemically rich molecular core (“hot corino”) could be detected at a distance of 4.9 kpc. Finally, the luminosity of the source (3000 $L_{\odot}$) is 3 orders of magnitude greater than that of a solar-type star.

Most of these problems can be solved if the disk is sufficiently inclined with respect to the line of sight. A luminosity of 3000 $L_{\odot}$ star corresponds to a stellar mass $M_{\ast}\simeq 9$ $M_{\odot}$ (Mottram et al. mott11 (2011)), which in turn implies an inclination \theta=\arcsin(\sqrt{\mbox{M_{\rm dyn}}/M_{\ast}})\simeq 19\degr. An even smaller $\theta$ is needed to justify $M_{\ast}>24$ $M_{\odot}$. While an almost face on disk is consistent with the findings of Molinari et al. (mol98 (1998)), the SE–NW orientation of the disk (and associated outflow) axis is not. In fact, the maps in Fig. 4 of Molinari et al. (mol98 (1998)), as well as the large-scale maps of Wu et al. (wu05 (2005)), suggest that the SiO flow is oriented NE–SW. In addition, the luminosity of a $>$24 $M_{\odot}$ star should exceed the observed luminosity of 3000 $L_{\odot}$. We thus believe that the model consisting of Keplerian disk with expansion is not acceptable.

4.2.2 Rotating disk

An alternative scenario is that the velocity pattern in Fig. 10 is due to the combination of a rotating disk oriented SE–NW, plus a bipolar outflow in the NE–SW direction. The SE–NW velocity gradient can be investigated through the position–velocity (PV) diagram along a suitable cut (dashed line in Fig. 10). Figure 12 shows this diagram for the lines of three species and in all cases the plot is vaguely reminiscent of the typical “butterfly-shaped” pattern consistent with Keplerian rotation (see the yellow curves in the figure; e.g. Cesaroni et al. cesa05 (2005)).

Under the assumption of Keplerian rotation, the velocity and position of the two putative “spurs” indicated in Fig. 12 can be used to estimate a dynamical mass, which for a Keplerian disk is the mass of the central object. We stress that we basically use the same approach as in Sect. 4.2.1 to derive an estimate of the stellar mass, with the only difference that here the disk is oriented SE–NW (and hence responsible for the corresponding velocity gradient), whereas in Sect. 4.2.1 was oriented NE–SW.

For a rotation velocity of $\sim$1 km s*-1* at a radius of $\sim$1$\aas@@fstack{\prime\prime}$5 or 0.036 pc, one obtains a stellar mass of $\sim$8 $M_{\odot}$ (see the Keplerian pattern in Fig. 12). Although this value is in good agreement with the core luminosity of 3000 $L_{\odot}$ estimated by Molinari et al. (2008b ), it is small compared to the mass of the core (24 $M_{\odot}$), inconsistent with the assumption of Keplerian rotation, which requires the gas mass to be negligible with respect to the stellar mass. This problem can be solved if the disk is sufficiently inclined. Since M_{\ast}=8\,\mbox{M_{\odot}}/\sin^{2}\theta (with $\theta$ angle between the disk axis and the line of sight), the stellar mass can exceed 24 $M_{\odot}$ for $\theta<35\degr$. However, the luminosity of a $>$24 $M_{\odot}$ star is $>7\times 10^{4}$ $L_{\odot}$ (Mottram et al. mott11 (2011)), much greater than the estimated value of 3000 $L_{\odot}$.

All the above assumes that the disk is Keplerian. If one drops this assumption and considers a self-gravitating disk, the dynamical mass estimated above (8 $M_{\odot}$) corresponds to the mass of the disk. As shown above, for a suitable inclination this can match the mass of 24 $M_{\odot}$ obtained from the millimeter continuum. Since the disk is self-gravitating, the stellar mass must be significantly less than the disk mass, consistent with $\sim$9 $M_{\odot}$ derived from the luminosity of 3000 $M_{\odot}$. In fact, a self-gravitating disk oriented close to face-on is in agreement with the findings of Molinari et al. (mol98 (1998)). Also, the size appears more consistent with those of “toroids” than to the typical diameters ($\sim$1000 au) of accretion disks around B-type stars (see Beltrán & de Wit beldew (2016)). Finally, a butterfly shaped PV diagram is not in contradiction with a self-gravitating disk because its rotation curve could mimic that of a Keplerian disk (e.g. Bertin & Lodato belo99 (1999), Douglas et al. doug13 (2013), Ilee et al. ile18 (2018)).

It is worth stressing that the uncertainties on the estimate of the dynamical mass are large. For example, assuming a radius of 2″ and a corresponding rotation velocity of 1.5 km s*-1*, one obtains 25 $M_{\odot}$. The core mass is also uncertain, as it is inversely proportional to the dust absorption coefficient (assumed 0.9 cm2g*-1* by BEU18), which might be underestimated by a factor $\sim$2 (see Table 1 of Ossenkopf & Henning osshen (1994)). Furthermore, the gas-to-dust mass ratio could be 100 instead of 150 (assumed by BEU18 in their calculations), which may cause an overestimate of the mass by a factor 1.5.

Despite all these caveats, we believe that the hypothesis of a self-gravitating disk oriented SE–NW and rotating about a $\sim$9 $M_{\odot}$ star is plausible. We remark that in this case, the NE–SW orientation of the disk axis would be consistent with that of the large-scale flow imaged by Wu et al. (wu05 (2005)).

4.2.3 Bipolar outflow

A third possibility to explain the observed velocity gradient is that of a compact bipolar outflow. Evidence for an outflow was provided also by Molinari et al. (mol98 (1998)) in terms of broad wings of the HCO+(1–0) and SiO(2–1) lines. Since no clear bipolar structure was identified, these authors proposed that the outflow is directed close to the line of sight. This is at odds with the small velocity difference ($\sim$1 km s*-1*) between the red- and the blue-shifted emission observed by us, because typical flow speeds are much greater ($>$10 km s*-1*) and projection effects should not matter if the flow is directed close to the line of sight. If an outflow is traced by the SiO emission, we believe that it should be very compact and with a significant inclination with respect to the line of sight.

To help distinguish between the disk and outflow scenarios we consider Fig. 9. Interestingly, the temperature rises from the center to the border of the core along the direction of the SE–NW velocity gradient, while the reverse occurs for the column density. This behaviour seems more compatible with a bipolar outflow than a rotating disk, because in the latter the temperature is expected to decrease with increasing radius. As a matter of fact, it would not be surprising if CH3CN was tracing an outflow, as several studies have detected CH3CN also in outflows (e.g. Leurini et al. leu11 (2011), Busquet et al. busq14 (2014), Palau et al. palau17 (2017)).

We can derive the parameters of the outflow, assuming that all the material traced by the CH3CN line participates in the expansion. As a first step, we estimate the abundance of CH3CN, $X_{\rm CH_{3}CN}$, from the ratio between the total number of CH3CN molecules and the total number of H2 molecules. The latter is obtained from the mass of 24 $M_{\odot}$ from Table 2, while the former can be computed by integrating the CH3CN column density over the map in Fig. 9. We obtain $X_{\rm CH_{3}CN}\simeq 10^{-10}$, an order of magnitude less than the typical value of hot molecular cores, but plausible for an object in an earlier evolutionary phase (see Gerner et al. gern14 (2014)). Then we calculate the momentum of the flow from the expression

[TABLE]

where $i$ is the inclination angle of the outflow axis with respect to the line of sight, the sum is extended over all the pixels of the map in Fig. 9, $N_{i}$ and $V_{i}$ are the CH3CN column density and velocity at pixel $j$, $V_{\rm sys}$ the systemic velocity ($\sim$ –50.7 km s*-1*), $\mu=2.8$ the mean molecular weight, $m_{\rm H}$ the mass of the hydrogen atom, $d$ the distance (4.9 kpc), and $\delta\Omega$ the solid angle of the pixel. We obtain $P\,\cos i\simeq 9$ $M_{\odot}$ km s*-1*. The corresponding timescale of the flow, $t_{\rm out}$, is given by the ratio between the SE–NW size of the region mapped in Fig. 9 and the maximum velocity range observed over the same region: t_{\rm out}\,\tan i=1\aas@@fstack{\prime\prime}9/1.8\leavevmode\nobreak\ \mbox{km\leavevmode\nobreak\ s{}^{-1}}=0.045\leavevmode\nobreak\ {\rm pc}/1.8\leavevmode\nobreak\ \mbox{km\leavevmode\nobreak\ s{}^{-1}}=2.4\times 10^{4} yr. From this one obtains the momentum rate $\dot{P}\,\cos^{2}i/\sin i=3.7\times 10^{-4}$ $M_{\odot}$ km s*-1* yr*-1*. Assuming momentum conservation in the flow, this value should be valid for the whole outflow, despite the fact that our estimate is obtained for the small region traced by the CH3CN emission.

Taken at face value and using the relationship between $\dot{P}$ and luminosity determined by Maud et al. (maud15 (2015)) for a distance-limited sample of outflows from massive YSOs (\mbox{\rm Log{}{10}}[\dot{P}]=-4.8+0.61\,\mbox{\rm Log{}{10}}[L]), the momentum rate implies a YSO luminosity of $\sim$175 $L_{\odot}$, far less than the estimate of $3\times 10^{3}$ $L_{\odot}$ obtained by Molinari et al. (2008b ), which instead corresponds to $\dot{P}=2.1\times 10^{-3}$ $M_{\odot}$ km s*-1* yr*-1*. To match this value, the outflow inclination must be $i\simeq 70\degr$. Such a large angle cannot be ruled out a priori, but is inconsistent with the claim of Molinari et al. (mol96 (1996)) that the outflow axis lies close to the line of sight. Moreover, the SE–NW orientation of the putative outflow in the plane of the sky is roughly perpendicular to that of the bipolar flow imaged by Wu et al. (wu05 (2005)) in the 12CO(1–0) line on the arcmin scale. While comparison between scales that differ by more than an order of magnitude must be taken with caution, we believe that the outflow interpretation for the observed CH3CN velocity field is probably the least likely of all those previously discussed by us.

Based on all the above, we conclude that the CH3CN emission in A1 is unlikely to trace an outflow and we are thus inclined to prefer the disk scenario. It is also possible that the observed velocity field is contributed by an accretion flow (through the disk), as suggested by the red-shifted self-absorption seen in the line profiles, as discussed later in Sect. 5.1.

4.3 Core A2

Figures 8 and 10 show that the LSR velocity increases slightly from A1 to A2. Figure 8 also shows that the velocity dispersion is significantly larger in A1 than in A2, while Figs. 4 and 6 demonstrate that the emission in all lines observed by us, as well as in the continuum, is weaker in A2. All of this supports our choice to distinguish this core from A1.

Another important feature denoting a difference between the two cores is the chemical richness. From Table 1, one sees that, while all listed species are detected in A1, only 2 of them are clearly revealed in A2. Finally, Fig. 9 demonstrates that both column density and temperature decrease going from A1 to A2, while the latter is half as massive as the former.

In conclusion, it seems unlikely that A2 is a starless core, given the relatively large temperature (typically $\sim$60 K) obtained from CH3CN. However, the lower mass and turbulence suggest that any YSO in A2 could be less massive and/or in an earlier evolutionary phase than in A1.

4.4 Core B

This core is clearly separated from A1 only in the continuum and CH3OH emission. Although B is a fainter emitter in most tracers with respect to A1 and A2, the mass and temperature are comparable to those of A1 and A2, hinting at a substantial similarity between these cores. The most striking feature of B is the ratio $M_{\rm H_{2}}/M_{\rm vir}$, which is the only one above unity. We will comment on this fact in more detail in Sect. 5.1, but we can already draw the conclusion that B must be virialized or on the edge of collapse and in all likelihood, this core is bound to evolve into a status analogous to that of A1 and A2. Indeed, the lack of emission in most of the observed lines (see Table 1) appears consistent with the idea that this core is in an early evolutionary phase.

4.5 Core C

The most interesting feature of core C is its possible association with the bipolar outflow described in Sect. 3.2. In Fig. 13 we show the PV diagram of the emission along the dashed line in Fig. 6, in three different lines. Several comments are in order for this plot. First of all, there is a velocity trend, outlined by the dashed line in the figure. Then it is worth noting that the H2CCO emission tracing core C is centrally located with respect to the SE and NW lobes (detected in CH3OH and H2CO) not only in space but also in velocity. Finally, the CH3OH and H2CO emission at the center is highly asymmetric in velocity, with a broad blue wing, evident also from the spectrum in Fig. 14. At the same time, the lack of an equally broad red wing could be explained by the presence of red-shifted self-absorption in a core undergoing infall. This is consistent with the blue-shifted “spot” visible in Fig. 8a, an indication of infall (see e.g. Mayen-Gijon et al. magi14 (2014)). These facts suggest that indeed the YSO powering the flow could be embedded in an infalling core characterized by large velocity dispersion. Whether this core is indeed core C (i.e. the one identified through the H2CCO emission) is unclear, as the latter is slightly offset from both the outflow axis (see Fig. 6) and the blue-shifted “spot”. In any case, we believe that this outflow cannot be associated with a very massive object, due to the weakness of the line emission from the lobes and the lack of a well defined core/YSO powering it.

4.6 Cores D and E

As one can see in Fig. 4, core D is only well defined in the CH3OH map, while other tracers present, at most, a tail of emission towards it. Similar to core D, also core E is clearly identified in only one tracer, this time the continuum emission. In summary, both cores have small masses and are recognizable as such only in one tracer. However, there is an important difference between the two: the velocity dispersion in D is much greater than that in E, as witnessed by the FWHM of the CH3OH line, which is 1.8 times larger in D than in E (see Table 2). We note that the line FWHM in D (3.6 km s*-1*) is even greater than that in A1 (3.2 km s*-1*), where various signposts indicate star formation activity.

5 Discussion

5.1 A cluster of massive stars in the making?

Based on the analysis presented above, we consider the possibility that IRAS 23385+6053 consists of multiple star-forming cores in an early stage of their evolution. Molinari et al. (2008b ) have provided evidence that this object is indeed in a protostellar phase, prior to the formation of an Hii region. However, their analysis assumes the presence of a single YSO, which we associate with A1, because this core seems to be the most active in terms of star formation and is hence likely to dominate the luminosity and mass estimates used by Molinari et al. (2008b ) to establish the evolutionary phase of IRAS 23385+6053.

We have shown that A1 is surrounded by other cores, most of which could contain embedded YSOs undergoing accretion. The first question we are asking is whether core formation proceeded through thermal Jeans fragmentation. This can be checked by comparing the Jeans length of the overall region with the mean separation of the cores. The latter is $\Delta\simeq 1\aas@@fstack{\prime\prime}2=5900$ au, whereas the former can be computed from $\lambda_{\rm J}({\rm au})=1.16\times 10^{6}\sqrt{T/n_{\rm H_{2}}}$, assuming the mean temperature ($T\simeq 30$ K) from CH3OH and deriving the gas density ($n_{\rm H_{2}}\simeq 2.3\times 10^{6}$ cm*-3*) from the ratio between the total core mass ($\sim$70 $M_{\odot}$ – see Table 2) and the volume occupied by the cores. The diameter of this volume (assuming spherical symmetry) is equal to the size ($\sim$4″) of the region over which the cores are distributed on the plane of the sky. We obtain $\lambda_{\rm J}\simeq 4200$ au. This value is only indicative, because in the pristine cloud both the density and temperature were lower than those at the present time. In consideration of these uncertainties, we believe that $\lambda_{\rm J}$ is comparable to $\Delta$ and thermal Jeans fragmentation could be a viable mechanism for the formation of the cores, consistent with the results other studies of massive cores at similar spatial scales (Palau et al. palau15 (2015, 2018); Ohashi oha18 (2018); Li et al. li19 (2019)).

The next question is whether the cores are in equilibrium. To address this issue, we consider the ratios $M_{\rm H_{2}}/M_{\rm vir}$ in Table 2. These are $<$1 in all cases except for core B, which seems to indicate that only one core could be effectively forming stars. To establish if this is indeed the case, we have searched for evidence of infall onto the cores. In Fig. 14 we show a comparison between the optically thick H2CO($3_{03}$–$2_{02}$) line and the thinner CH3OH($4_{2}$–$3_{1}$) E1 line, towards all the cores. We stress that the H2CO spectra are obtained by merging the NOEMA and 30-m data to also recover the extended emission. Therefore, the dips in the profiles cannot be ascribed to missing flux filtered out by the interferometer.

The H2CO transition presents red-shifted self-absorption towards almost all cores, except B. In fact, the H2CO profiles are skewed towards blue-shifted velocities with respect to the CH3OH line, whose peak velocity we take as a proxy for the systemic velocity of the core. This indicates the presence of red-shifted self-absorption, a well-known indicator of infall (see e.g. Mardones et al. mardo (1997) and references therein), which in turn suggests that the cores could be forming stars in the main accretion phase. This result is apparently inconsistent with the fact that all cores except B have $M_{\rm H_{2}}<M_{\rm vir}$. We propose that, instead, the presence of infalling material could exert an external pressure on the cores, sufficient to stabilize them and even trigger their collapse. This scenario could explain why some of the cores are already forming stars, despite the low value of $M_{\rm H_{2}}/M_{\rm vir}$.

In practice, we can take into account the effect of external pressure using Eq. (8) of Field et al. (field (2011)), where the surface density is expressed as $M_{\rm H_{2}}/(\pi\,R_{\rm c}^{2})$. The virial mass is obtained from this equation as

[TABLE]

where $R_{\rm c}$ is the core radius, $\Delta V$ the line FWHM, $G$ the gravitational constant, $p_{\rm e}$ the external pressure, $p_{0}=\frac{15}{16\pi}\frac{(\Delta V)^{4}}{(8\ln 2)^{2}\,G\,R_{\rm c}^{2}}$, and we assume that the core has constant density (i.e. $\Gamma=3/5$ in the notation of Field et al. field (2011) – see also MacLaren et al. mclar (1988)). If $p_{\rm e}>p_{0}$ no equilibrium configuration can be attained and the core collapses.

The external pressure on the core is the sum of the thermal pressure of the envelope plus the ram pressure of the infalling gas444 Our approach of including the ram pressure in the external pressure on the core surface is not strictly correct in the hydrostatic scenario of Field et al. (field (2011)). However, it has been shown by Ballesteros-Paredes et al. (bapa11 (2011)) that such an approach yields results comparable to those obtained from the correct treatment of infall. , namely $p_{\rm e}=n_{\rm H_{2}}\,(k\,T+\mu\,m_{\rm H}\,\varv_{\rm inf}^{2})$, with $k$ Boltzmann constant, $T$ gas temperature, and $\varv_{\rm inf}$ infall velocity. We assume $\varv_{\rm inf}\simeq 2$ km s*-1* from the difference between the emission peak and absorption dip in the spectra of Fig. 14, while $n_{\rm H_{2}}\simeq 2.3\times 10^{6}$ cm*-3* and $T=30$ K are estimated as previously explained. From these and the parameters in Table 2, we find that $p_{\rm e}\simeq 9.7\times 10^{-7}$ dyn cm*-2* and $p_{0}\simeq 5.4\times 10^{-8}$–$9.4\times 10^{-7}$ dyn cm*-2*. Therefore, in all cases $p_{\rm e}>p_{0}$ and the cores appear to be bound to collapse.

Under the assumption that all cores are collapsing, we compute the accretion rates, $\dot{M}_{\rm acc}$, and corresponding accretion timescales, $t_{\rm acc}$, from the expressions $\dot{M}_{\rm acc}=[3/(8\,\ln 2)]^{3/2}\Delta V^{3}/G)$ and $t_{\rm acc}=M_{\rm H_{2}}/\dot{M}_{\rm acc}$, where $G$ is the gravitational constant. The values are listed in Table 2. While $\dot{M}_{\rm acc}$ is consistent with the values for stars $>$10 $M_{\odot}$ (see Beltrán & de Wit beldew (2016)), $t_{\rm acc}$ seems quite short compared to the typical timescale for the formation of a massive star (\sim$$10^{5} yr) and would imply that the cores are short lived or very young. However, on the one hand not all the core material is necessarily accreting onto a star, and on the other hand the cores can be loading fresh material from the surrounding envelope, as indicated by the evidence for infall discussed above. Our estimate of $t_{\rm acc}$ is hence to be considered a lower limit.

Another topic that is worth considering is the chemical richness of the cores. A detailed analysis of the core chemical content requires careful identification of the molecular lines detected in all cores. This is a complex process that goes beyond the scope of our study. However, we adopt a different approach, estimating the amount of molecular line emission in the different cores. For each core we consider the spectrum obtained with the broad-band correlator, WIDEX, and – following Cesaroni et al. (cesa17 (2017)) – we calculate the fraction of spectral channels, $f$, where emission is detected above 5$\sigma$. In Fig. 15, this parameter is plotted versus the corresponding mean S/N, obtained as the ratio between the mean flux over the detected channels and the 1$\sigma$ RMS of the spectrum. In this plot, for the same value of S/N, sources with higher values of $N_{\rm det}/N_{\rm tot}$ are more line rich and, presumably, also more chemically rich. Therefore, we can state that D is less rich than A1, and A2 is richer than C. As for B and E, they appear line poor, but their S/N is worse than for the others and it is thus possible that part of the line emission is not detected only because it is weaker than in the other cores.

With all the above in mind, we can conclude that some chemical diversity is present in the IRAS 23385+6053 region. However, the spread in $f$ among the cores is significantly less than that found by Cesaroni et al. (cesa17 (2017)) in their sample of hot molecular cores (see their Table 3). While in their case the values of $f$ are an order of magnitude greater, they span a wider range (a factor 4.6) than in our sample (a factor 2.1). This yields a twofold result: on the one hand, the chemistry in our cores is less evolved than in typical hot molecular cores; on the other hand, our cores do not differ much each other, which suggests that they might form the same type of stars.

In conclusion, we believe that our results provide convincing evidence that the multiple cores in IRAS 23385+6053 represent a massive cluster555 With “massive cluster” here we mean a stellar cluster with at least one massive (O-type or early B-type) YSO. in the making. Although quite similar in mass (apart from C), the cores appear to be in (slightly) different evolutionary stages, given their differences in terms of temperature and velocity dispersion.

Is the purported cluster gravitationally bound? A rough estimate of its stability can be obtained from the virial theorem. We assimilate the core cluster to a bunch of mass particles, characterized by the masses and line-of-sight velocities given in Table 2. From the standard deviation of the velocity distribution ($\sim$0.49 km s*-1*) and the maximum separation among the cores ($\sim$3″), one obtains an equivalent virial mass of $\sim$10 $M_{\odot}$, to be compared with the total mass of the cores of $\sim$70 $M_{\odot}$. Despite the large uncertainties of the method, the difference between the two masses appears large enough to indicate that the core cluster is in all likelihood gravitationally bound.

5.2 Comparison with other high-mass star-forming regions

The scenario depicted above is reminiscent of the case of W33A, investigated by Maud et al. (maud17 (2017)). Both sources, when observed at higher resolution, reveal a cluster of cores, possibly in different evolutionary phases, and do not present iron-clad evidence of Keplerian rotation, although the presence of disks on smaller scales is suggested in W33A by detailed modelling (Izquierdo et al. izq18 (2018)). Also, in W33A accretion onto the cores appears to proceed through a spiral-like filament, while in our case the line profiles hint at infall onto most of the cores. Despite these similarities, there are also differences between IRAS 23385+6053 and W33A, the most important of which are the spatial scales investigated in the two cases. The spatial resolution is $\sim$10 times worse in IRAS 23385+6053 ($\sim$4400 au) than in W33A ($\sim$500 au), and the cores are distributed over a region of $\sim$3500 au in W33A as opposed to $\sim$20000 au in IRAS 23385+6053. Moreover, W33A is $\sim$10 times more luminous but $\sim$20 times less massive, and seems more chemically rich than IRAS 23385+6053. All these features reinforce our conviction that IRAS 23385+6053 is a young massive cluster in the making, with a (proto)stellar content yet in an early evolutionary phase (as predicted by Molinari et al. mol98 (1998); 2008b ), whereas W33A represents a more evolved stage.

More similar to our source is probably IRAS 05358+3543, where at least 4 protostellar cores have been identified by Beuther et al. (beu07 (2007)). Although the spatial resolution is $\sim$900 au, better than in our case, the core cluster spans a region of $\sim$16000 au, alike to IRAS 23385+6053. The bolometric luminosity (\sim$$6\times 10^{3} $L_{\odot}$) and total mass of all cores ($\sim$34 $M_{\odot}$; see Table 3 of Beuther et al. beu07 (2007)) are also similar to those of our source. Despite these similarities, the two objects also present two substantial differences: no free-free continuum emission has been detected from the cores in IRAS 23385+6053, whereas a hypercompact Hii region is found in IRAS 05358+3543; and multiple outflows are known to be associated with the cores in this source, unlike our case where only questionable evidence of two bipolar outflows has been found. These differences seem to suggest that IRAS 05358+3543 is in a more advanced stage with respect to IRAS 23385+6053. However, one should keep in mind that no deep, sub-arcsecond continuum observation at centimeter wavelengths has been performed towards IRAS 23385+6053. It is worth noting that the presence of the two Hii regions shown in Fig. 1 could make it difficult to detect faint free-free continuum emission from any putative hypercompact Hii region deeply embedded in one of the cores.

6 Summary and conclusions

In the context of the CORE project (P.I. Henrik Beuther) we have performed observations of the continuum and molecular line emission at $\sim$1.4 mm from the high-mass star-forming region IRAS 23385+6053. This is believed to be the youngest of the whole CORE sample. We have revealed six cores and possibly an outflow from a low-mass star. One of the cores (A1) could contain a self-gravitating disk rotating about a $\sim$9 $M_{\odot}$ star, responsible for most of the 3000 $L_{\odot}$ estimated for IRAS 23385+6053. The other cores are less massive and colder, but all of them appear to be on the edge of collapse. The accretion rates estimated from the velocity dispersions are consistent with those typical of high-mass star-forming regions. We conclude that IRAS 23385+6053 contains a sample of cores that are forming, or are bound to form, a cluster of massive stars.

Acknowledgements.

It is a pleasure to thank Daniele Galli for stimulating discussions and Sergio Molinari for critically reading the manuscript. We thank the IRAM technical staff for their support in this project. DS acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: “Building a Habitable Earth” (SE 1962/6-1). HB, AA, JCM, and SS acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648505. RK 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. AP acknowledges financial support from UNAM-PAPIIT IN113119 grant, México. RGM acknowledges support from UNAM-PAPIIT Programme IN104319. This work is based on observations carried out under project number L14AB008 with the IRAM interferometer and 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). XCLASS development is supported by BMBF/Verbundforschung through the Projects ALMA-ARC 05A11PK3 and 05A14PK1 and through ESO Project 56787/14/60579/HNE.