Observations of high-order multiplicity in a high-mass stellar protocluster

Abstract

The dominant mechanism forming multiple stellar systems in the high-mass regime (M_* ≳ 8 M_⊙) remained unknown because direct imaging of multiple protostellar systems at early phases of high-mass star formation is very challenging. High-mass stars are expected to form in clustered environments containing binaries and higher-order multiplicity systems. So far only a few high-mass protobinary systems, and no definitive higher-order multiples, have been detected. Here we report the discovery of one quintuple, one quadruple, one triple and four binary protostellar systems simultaneously forming in a single high-mass protocluster, G333.23–0.06, using Atacama Large Millimeter/submillimeter Array high-resolution observations. We present a new example of a group of gravitationally bound binary and higher-order multiples during their early formation phases in a protocluster. This provides the clearest direct measurement of the initial configuration of primordial high-order multiple systems, with implications for the in situ multiplicity and its origin. We find that the binary and higher-order multiple systems, and their parent cores, show no obvious sign of disk-like kinematic structure. We conclude that the observed fragmentation into binary and higher-order multiple systems can be explained by core fragmentation, indicating its crucial role in establishing the multiplicity during high-mass star cluster formation.

Main

High-mass stars in the Milky Way are overwhelmingly (>80%; refs. ^1,2,3,4,5,6) found in binaries or higher-order multiplicity systems that play a key role in governing cluster dynamics and stellar evolution^7,8. However, it is yet unclear whether they are predominantly formed from in situ fragmentation at various scales (for example, disks^9,10, cores^11,12 or filaments¹³) or subsequent stellar capture in clusters¹⁴ because direct measurements of their initial configuration and properties at the early phases of cluster formation have been unattainable^{15,16,17,18,19,20,21,22,23,24,25,26,27}.

We report the direct imaging of one quintuple, one quadruple, one triple and four binary systems in the high-mass protocluster G333.23–0.06 (hereafter G333) using Atacama Large Millimeter/submillimeter Array (ALMA) long-baseline observations (Fig. 1 and Extended Data Fig. 1). These binary and higher-order systems are detected in the high-resolution (angular resolution of θ ≈ 0.05″, equivalent to 260 AU at the source distance of 5.2 kpc; ref. ²⁸) 1.3 mm dust continuum image. The detected condensations have radii between 153 and 678 AU (Table 1). We refer to both binary and higher-order systems simply as multiple systems in what follows, and we only make an explicit distinction between the two when necessary.

Fig. 1: Continuum images of ATCA (3.3 mm), ALMA low-resolution (1.3 mm) and ALMA high-resolution (1.3 mm) observations.

a, ATCA 3.3 mm continuum image (θ ≈ 2.42″). The contour levels are [5, 8, 11, 14, 17, 20, 23, 26, 29] × σ, where root mean squared (rms) noise of continuum image is σ = 0.6 mJy per beam. b, ALMA low-resolution (θ ≈ 0.32″) 1.3 mm continuum image. The contour is the 7σ, where σ = 0.16 mJy per beam. The ellipses show the identified cores based on the low-resolution continuum image. The red pluses are the Class II CH₃OH maser⁴². c–g, ALMA high-resolution (θ ≈ 0.05″) 1.3 mm continuum image for multiple systems. The green crosses present the identified condensations based on the ALMA high-resolution continuum image. The contour is the 7σ, where σ = 0.05 mJy per beam. The dense cores (no. 1, no. 2, no. 3, …) and condensations (C1, C2, C3, …) are numbered in order of descending integrated intensity. c, Dense core no. 2 fragments into quadruple condensation system (C8, C10, C14, C17) and dense core no. 3 fragments into triple condensation system (C6, C12, C26). d, Dense core no. 1 fragments into quintuple condensation system (C1, C3, C4, C5, C16) and dense core no. 7 fragments into binary condensation system (C11, C29). e, Dense core no. 17 fragments into binary condensation (C22, C38). f, Dense core no. 30 fragments into binary condensation (C39, C42). g, Dense core no. 15 fragments into binary condensation (C35, C40). The white ellipses in the lower left corner of each panel denote the synthesized beam of continuum images. ICRS, International Celestial Reference System.

Table 1 Properties of condensations

The projected separations of these multiple systems are between 327 and 1,406 AU, with a mean value of 731 AU (Extended Data Fig. 2), in good agreement with the typical projected separation of 700 AU in the simulation of multiple star formation via core fragmentation²⁹. The ambient gas masses (M_amb) of these multiple systems range from 0.19 to 1.47 M_⊙ on the basis of the thermal dust emission (Methods). These masses are regarded as lower limits because the observations suffered from missing flux in the interferometer data. We focus primarily on the multiple systems that are embedded in a single dense core (typical radius of ∼2,100 AU; Fig. 1). The quintuple system consists of a small group of condensations (C1–C4–C5–C16), which is tightly connected as seen in dust continuum emission, and a condensation (C3) slightly separated. That is nevertheless part of the original parental core (Fig. 1). The quadruple system includes two binary configurations (C10–C14 and C8–C17), and the triple system composes three slightly separated condensations (C6, C12 and C26). The binary systems are C11–C29, C22–C38, C39–C42 and C35–C40.

In addition to the projected proximity on the sky of observed condensations, the line-of-sight velocity is another important diagnostic tool to determine whether members of a multiple system are physically associated. All members of each multiple system have similar centroid velocities (Extended Data Fig. 3 and Methods), indicating that the members are physically associated and share a common origin.

Multiple systems formed via core fragmentation

The parent cores of the multiple systems are revealed in the lower-resolution (θ ≈ 0.32″, equivalent to 1,664 AU) 1.3 mm dust continuum image (Fig. 1 and Extended Data Fig. 4). The dense cores are identified with radii ranging from 927 to 3,443 AU (Table 1). The multiple condensation systems are embedded in the dense cores. Each condensation likely harbours embedded protostellar object(s) as evidenced by the presence of hot (T_gas = 108–532 K) and warm gas resulting from internal heating (Methods), except for C39 and C42 where no significant molecular warm transitions (that is, upper energy level of E_u/k > 45 K) are detected.

There is no obvious sign of a disk-like kinematic structure around any of the multiple systems and their parent cores in any of the lines we examined (Methods), including the typical disk tracers, for example, CH₃CN, ¹³CH₃CN and CH₃OCHO (Fig. 2 and Extended Data Figs. 5–7), and dense gas tracers, for example, CH₃OH, ¹³CH₃OH, SO₂, SO, HC₃N, HNCO, NH₂CHO, H₂CO and H${}_{2}^{13}$CO. Meanwhile, the presence of the SiO outflows indicates that the accretion disk is not viewed face-on, which rules out the scenario of a weak velocity gradient resulting from projection of a face-on geometry (Fig. 2 and Extended Data Fig. 8). In addition, the scenario in which multiple systems form by dynamical capture in a forming cluster, which have typical separations >10³ AU and significant different velocities³⁰, is inconsistent with the observed separation distributions and small velocity differences.

Fig. 2: Examples of intensity-weighted velocity maps and position-velocity (PV) diagrams of molecular lines for the quintuple system.

Left, middle and right columns present the results derived from the CH₃CN 12₄ − 11₄, ¹³CH₃CN 13₃ − 12₃ and CH₃OCHO 20_0,20 − 19_0,19, respectively. a–c, We present the intensity-weighted velocity maps derived from the ALMA high-resolution data of CH₃CN 12₄ − 11₄ (a), ¹³CH₃CN 13₃ − 12₃ (b), and CH₃OCHO 20_0,20 − 19_0,19 (c). The black and magenta ellipses show the condensations and their parent cores, respectively. The red plus marks the Class II CH₃OH maser position. The blue and red arrows show the directions of the main outflow seen in the SiO emission from the ALMA low-resolution data. The grey contours show the velocity-integrated intensity maps with levels of [5, 10, 15, 20] × σ_rms, where rms noise of velocity-integrated intensity map is σ_rms = 7.5 mJy beam⁻¹ km s⁻¹. The black ellipses in the lower left corner of each panel denote the synthesized beam of continuum images. The remaining velocity maps are presented in Extended Data Figs. 5–7. d–f, We show the PV diagram maps derived from the high-resolution data of CH₃CN 12₄ − 11₄ (d), ¹³CH₃CN 13₃ − 12₃ (e) and CH₃OCHO 20_0,20 − 19_0,19 (f). Contour levels start at 4σ_rms and increase in step of 1σ_rms interval, where σ_rms is 1.2 mJy per beam. The cut lines of the PV diagram are indicated in a–c with black lines. The blue dashed lines in the vertical and horizontal directions show the position of dense core no. 1 and its systemic velocity of V_lsr = –85 km s⁻¹.

If the parental disks are small and the separations of fragments are significantly widened within a short timescale³¹, we cannot completely rule out the possibility that disk fragmentation is responsible for the close binary systems (for example, C22–C38 and C35–C40). Overall, these results demonstrate that the majority of detected multiple systems is formed from core fragmentation, although disk fragmentation may still occur on smaller scales than those we can resolve with the current spatial resolution. The measured separations of the multiple systems (a mean value of 731 AU) are smaller than the expected Jeans lengths of 5,000–10,000 AU for thermal Jeans fragmentation of the parent cores of the multiple systems (Methods). The estimated Jeans lengths are conservative upper limits to the separations of fragments since the derived volume densities of parental cores are lower limits due to the missing flux. If one assumes that the turbulence is acting as an isotropic support, the turbulent Jeans fragmentation yields even larger separation than that of the thermal one because the non-thermal velocity dispersion is higher than the thermal velocity dispersion. On the other hand, the anisotropic turbulent motions could promote collapse³², in which the Jeans length could be reduced. The discrepancy indicates that the core fragmentation, which facilitates the formation of multiple system, is not merely regulated by the thermal pressure and/or turbulence, but additional mechanisms might also be important. For example, ongoing global collapse and dynamical interactions of multiple systems could lead to inward migration³³, which moves the fragments closer together.

Masses of the central protostars

Conventionally, the mass (M_*) of the central protostar can be estimated through modelling the rotation in a Keplerian disk. However, the estimation of dynamical masses of the central protostars is prevented by the non-detection of disk kinematic structures toward the multiple protostellar systems. On the other hand, a roughly M_* can be estimated from the bolometric luminosity and a zero-age main sequence (ZAMS) assumption for the young protostars, which could provide a mass comparable to the dynamical mass obtained from Keplerian rotation³⁴. The luminosities of the central protostars are estimated from the temperature profile with the assumption that the spectral index of dust emissivity at far-infrared wavelengths is β = 1 (Methods). The derived luminosities range between 50.9 and 1.5 × 10⁵ L_⊙, corresponding to B9–O6.5 spectral type ZAMS stars³⁵. The masses are estimated to be between roughly 2.5 and 26.1 M_⊙ according to the mass–luminosity relation³⁶ (Table 1). The derived masses are regarded as upper limits because the luminosities could be overestimated (Methods). Three binary systems (C22–C38, C39–C42, C34–C40), C12 and C29 do not have temperature measurements due to non-detection of CH₃CN and its isotopologues. There are at least three condensations (C1, C10, C14) with a mass of >8 M_⊙, even if the luminosity reaches a minimum value when β decreases down to 0. This shows that indeed high-mass protostars exist in the quintuple and quadruple systems, which is consistent with the presence of Class II CH₃OH maser emission toward these regions (Fig. 1). The results indicate that high- and low-mass multiple protostellar systems are simultaneously forming within G333.

Stability of the multiple systems

To determine whether a multiple system is bound, we used the ambient mass M_amb to compute the kinetic (E_i) and gravitational (W_i) energies for each member of the multiple systems (Methods). The derived E_i is smaller than |W_i| for all condensations (Fig. 3), suggesting that these multiple systems are gravitationally bound, except for two condensations (C10 and C14). If including the central protostellar masses, the multiple systems would be even more gravitationally bound because the gravitational energy will be even higher. Indeed, the E_i and W_i computed from protostellar mass M_* show that the E_i/|W_i| ratio is below 0.1 and smaller than those derived from the ambient mass (Fig. 3). In addition, the E_i/|W_i| ratio will be even smaller when both M_amb and M_* are included. Therefore, the multiple systems are gravitationally bound at the present stage. With 20 identified multiple systems in G333, the observed multiplicity fraction is MF = 20/44 ≈ 45%, which is the fraction of systems that are multiples (binary, triple and so on), and companion frequency is CF = 46/44 ≈ 1.0, which is the average number of companions per system. The derived MF and CF are higher than those measured in Orion and Perseus star-forming regions for a similar separation range of 300–1,400 AU (ref. ³⁷), indicating that the multiplicity could be higher in denser cluster-forming environments. The estimated MF and CF are regarded as lower limits because further fragmentation might occur at smaller scales than what we can resolve with the current observations, and low-mass objects could be missed due to the limited sensitivity. The results indicate that the multiplicity in clusters is established in the protostellar phase.

Fig. 3: The kinetic-to-gravitational energy ratio E_i/|W_i| as a function of mass for the multiple systems.

The multiple systems with kinetic-to-gravitational energy ratio below unity are considered to be gravitationally bound. The circles and stars symbols are the results derived from ambient mass (M_amb) and protostellar mass (M_*), respectively. E_i/|W_i| has been estimated with four different methods: (1) line-of-sight velocity difference and on-sky separation (refer to one-dimensional, 1D; black symbols), (2) three-dimensional (3D) velocity difference ($\sqrt{3}$ times the line-of-sight velocity difference) and on-sky separation (red symbols), (3) line-of-sight velocity difference and 3D separation ($\sqrt{2}$ times the on-sky separation; blue symbols) and (4) 3D velocity difference and 3D separation (orange symbols). The black dashed line marks E_i/|W_i| = 1. The members of the quintuple system are marked with different colour shadows, that is, C1 (blue), C3 (cyan), C4 (red), C5 (yellow) and C16 (black). There are two condensations (C10 and C14) with E_i/|W_i| > 1 for the 3D velocity scenario in the case of using M_amb. If the central protostar mass is considered, the E_i/|W_i| of these two condensations is smaller than 1. This figure indicates that all multiple systems are gravitationally bound.

Perspectives

The discovery of these quintuple, quadruple, triple and binary protostellar systems is the best observational evidence to show the imprints of core fragmentation in building multiplicity in high-mass cluster-forming environments. Although we cannot test if disk fragmentation is more important at smaller scales than what we have observed so far, we expect that more systems similar to G333 will be discovered given the high resolutions and high sensitivities of ALMA observations. The statistics of these systems will help to benchmark the relative contribution of core fragmentation to the population of multiple stars in high-mass star clusters. Their properties will determine the initial conditions of multiple system formation, as well as the dynamical evolution in a cluster environment.

Methods

High-mass star-forming region G333.23–0.06

G333.23–0.06 is a typical high-mass star-forming region^38,39,40,41 at a distance of 5.2 kpc (ref. ²⁸) associated with Class II CH₃OH maser emission⁴², which can only be excited in high-density regions by strong radiation fields, making it exclusively found in high-mass star-forming regions^43,44. G333 has a mass reservoir of ∼3,000 M_⊙ with a mean column density of 1.6 × 10²³ cm⁻² within a 1.2 pc radius⁴⁵. As a representative example of high-mass star-forming region, G333 provides an ideal laboratory to study the binary and high-order multiple system formation in the cluster environment.

Observations and data reduction

Observations of G333 were performed with ALMA in Band 6 (at the wavelength of 1.3 mm) with the 12 m array using 41 antennas in configuration similar to C40-5 (hereafter short-baseline or low-resolution) on 05 November 2016 and 42 antennas in configuration similar to C43-8 (hereafter long-baseline) on 28 July 2019 (Project ID: 2016.1.01036.S; PI: Sanhueza). Observations were obtained as part of the Digging into the Interior of Hot Cores with ALMA (DIHCA) project^27,46,47. The baseline lengths are 18.6–1,100 m and 91–8,547 m for short- and long-baseline observations, respectively. The correlators were tuned to cover four spectral windows with a spectral resolution of 976.6 KHz (∼1.3 km s⁻¹) and a bandwidth of 1.875 GHz. These windows covered the frequency ranges of 233.5–235.5 GHz, 231.0–233.0 GHz, 219.0–221.0 GHz and 216.9–218.7 GHz. The quasar J1427–4206 was used for flux and bandpass calibration, and J1603–4904 for phase calibration. The total on-source time toward the G333 is 6 minutes for short-baseline observations and 19.6 minutes for long-baseline observations. The phase centre used is (α(ICRS), δ(ICRS)) = 16 h 19 m 51.20 s, −50°15′13.00″.

The visibility data calibration was performed using the CASA (v.5.4.0-70) software package⁴⁸. We produced continuum data from line-free channels and continuum-subtracted data cubes for each observation epoch using the procedure described in ref. ⁴⁶. We performed phase-only self-calibration using the continuum data and the self-calibration solutions were applied to data cubes. To recover the extended emission, we combined the short-baseline and long-baseline self-calibration data for both continuum and data cubes (hereafter combined or high-resolution data). We produced images for short-baseline and combined data sets, separately. We used the TCLEAN task with Briggs weighting and a robust parameter of 0.5 to image the continuum. The resultant continuum images have a synthesized beam of 0.35″ × 0.30″ (1,820 AU × 1,560 AU, Fig. 1b) with a position angle of PA = −46.18°, 0.059″ × 0.038″ (307 AU × 198 AU) with a PA = 56.23° and 0.066″ × 0.039″ (343 AU × 203 AU, Fig. 1c–g) with a PA = 54.47° for short-baseline, long-baseline and combined dataset, respectively. The achieved 1σ rms noise levels continuum images are about 0.16 mJy per beam, 0.05 mJy per beam and 0.05 mJy per beam for short-baseline, long-baseline and combined data, respectively. Data cubes for each spectral window were produced using the automatic masking procedure YCLEAN⁴⁹, which automatically cleans each map channel with custom-made masks. More details on the YCLEAN algorithm can be found in ref. ⁴⁹. The lines images 1σ rms noise are about 10 mJy per beam, 3 mJy per beam and 3 mJy per beam with a channel width of ∼0.65 km s⁻¹ for short-baseline, long-baseline and combined data, respectively. The largest recoverable angular scales are 3.5″ for short-baseline and combined data, as determined by the short baselines in the array.

The Australia Telescope Compact Array (ATCA) 3.3 mm continuum image is retrieved from ref. ⁴¹ (Fig. 1a). All images shown in this article are before primary beam correction, while all measured fluxes have the primary beam correction applied.

Dense core and condensation identification

To describe the dense molecular structures, we follow the nomenclature in the literature in which cores refer to structures with sizes of ∼10³−10⁴ AU, and condensations refer to substructures within a core.

We use the astrodendro algorithm and CASA-imfit task to extract dense cores from short-baseline 1.3 mm continuum image and condensations from the combined and long-baseline only 1.3 mm continuum images. The astrodendro identifies the changing topology of the surfaces as a function of contour levels and extracts a series of hierarchical structure over a range of spatial scales⁵⁰. The performance of astrodendro in characterizing the dense structure parameters (for example, size and position angle) is not always good, while CASA-imfit performs better in this regard via a two-dimensional Gaussian fit to the emission.

Therefore, we used astrodendro to preselect dense structures (that is, the leaves in the terminology of astrodendro) from the 1.3 mm continuum images. We then use the parameters of the preselect structures from astrodendro as input to CASA-imfit for more accurate measurement of their parameters, including peak position, peak flux (I_peak), integrated flux (S_ν), major and minor axis sizes (full-width at half-maximum; FWHM_maj and FWHM_min) and position angle (PA).

The following parameters are used in computing the dendrogram: the minimum pixel value min₋value = 5σ, where σ is the rms noise of the continuum image; the minimum difference in the peak intensity between neighbouring compact structures min₋delta = 1σ; and the minimum number of pixels required for a structure to be considered an independent entity min₋npix = N, where N is the number of pixels in the synthesized beam area.

To remove suspicious condensations around the strong emission regimes caused by the diffuse emission in the combined data, we have performed a cross-comparison of condensation catalogue derived from the combined data with the condensations revealed by the long-baseline only data. We identified 30 dense cores in the short-baseline 1.3 mm continuum image and 44 condensations in the combined 1.3 mm continuum image. Extended Data Figs. 4 and 1 show the identified dense cores and condensations, respectively (see also Table 1 and Extended Data Fig. 2 for the properties of multiple systems).

Centroid velocity of condensation

The centroid velocity (V_lsr) of each condensation is determined by Gaussian fitting a CH₃OH 4_2,2 − 3_1,2 (E_u/k = 45.46 K) line that is detected in the majority of condensations to measure V_lsr in the same manner. The measured V_lsr have been validated by comparing with other dense gas tracers. We identify no clear velocity difference between the members of each multiple system (Table 1), that is, ΔV_lsr < 1 km s⁻¹ that is smaller than the line-of-sight velocity differences (2.0–9.5 km s⁻¹) of the binary protostars in refs. ^22,24, indicating that all members of each multiple system are associated with the same region. We show the moment maps of CH₃OH in Extended Data Fig. 3 for reference.

Jeans length

If the fragmentation is governed by Jeans instability, the Jeans length λ_J, which is the separation between fragments, can be calculated by⁵¹

$${\lambda }_{J}=\sqrt{\frac{\uppi \,{\sigma }_{{{{\rm{eff}}}}}^{2}}{G\,\rho }}=0.06{\,{\rm{pc}}}{\times}\left(\frac{{\sigma }_{{{{\rm{eff}}}}}}{0.188\ {{{\rm{km}}\,{\rm{s}}^{-1}}}}\right){\left(\frac{\rho }{1{0}^{5}\ {{{{\rm{cm}}}}}^{-3}}\right)},$$

(1)

where σ_eff is the effective velocity dispersion, ρ is the volume density and G is the gravitational constant. The σ_eff equals the sound speed c_s for thermal Jeans fragmentation. The temperatures and volume densities of the parent cores used in Jeans analysis are T = 80–340 K and ρ = 3 × 10⁶–1 × 10⁷ cm⁻³, respectively. The derived thermal Jeans lengths range from 5,000 to 10,000 AU. In the turbulent Jeans fragmentation scenario, the σ_eff includes thermal and non-thermal velocity components, ${\sigma }_{{{{\rm{eff}}}}}=\sqrt{{c}_{{{{\rm{s}}}}}^{2}+{\sigma }_{{{{\rm{nth}}}}}^{2}}$, under the assumption that the turbulence is acting as an isotropic of support. The measured line widths of CH₃CN are 3–5 km s⁻¹ for the parent cores, resulting turbulent Jeans lengths of 20,000–54,000 AU.

Search for disk kinematic structure

The observations cover the typical disk tracers, including CH₃CN and its isotopologues, and CH₃OCHO, as well as other dense gas tracers, for instance H₂CO and its isotopologues, CH₃OH and its isotopologues, HC₃N, NH₂CHO, SO₂, SO, HNCO, HCOOH, ¹³CS and OCS. Using these molecular lines, we have searched for disk-like rotating structures for the multiple systems and their parent cores with both short-baseline and combined data that have a channel width of ∼0.65 km s⁻¹. The dense gas tracers are not sufficiently strong to allow a reliable determination of kinematic information for three binary systems (C22–C38, C39–C42, C35–C40) in the combined data.

There is no obvious sign of disk kinematic structures toward the parent cores of multiple systems in any of the lines we examined based on short-baseline and combined data (Fig. 2 and Extended Data Fig. 5). There are some lines with a velocity gradient in some dense cores, but no clear Keplerian disk-like rotating structures are found in the position-velocity (PV) diagram toward these cores. The velocity gradients trace either the outflows or the large-scale gas motions (for example, gas flow, toroidal motions⁵²). These dense cores are associated with unipolar, bipolar and/or perpendicular outflows identified by the SiO emission from the ALMA short-baseline data (Extended Data Figs. 5 and 8). The detected misaligned outflows indicate that the embedded multiple systems do not come from the same co-rotating structures⁵³. This further suggests that the quintuple, quadruple and triple systems are formed from core fragmentation⁵⁴. The detailed analysis of molecular outflows is beyond the scope of this article and will be presented in a future paper.

We examined the multiple systems following the same routine but using the combined data and similarly found signs of velocity gradient in some condensations, but no obvious rotational signatures of disks (Fig. 2 and Extended Data Fig. 6). Some velocity gradients are likely dominated by the outflows, while the others require higher angular resolution and sensitivity to spatially resolve the origin (for example, unresolved outflows or accretion flows). At the early evolutionary stages of massive star formation, the size of disks might still be very small, potentially smaller than 100 AU (refs. ^55,56,57). Considering the G333 is still at early evolutionary phases, the non-detection of disk structures indicates that the size and/or mass of disks are smaller than what we can resolve with the current observations. The current spatial resolution is ∼260 AU and 3σ point source mass sensitivity of ∼0.03 M_⊙ assuming a temperature of 50 K.

As shown in Extended Data Fig. 6, there is a redshifted velocity feature surrounding the blueshifted velocity toward C10 and C14. Two velocity components are detected toward C10 and C14. We have inspected these two velocity components separately and found no obvious disk kinematic (Extended Data Fig. 7). Several mechanisms could lead to two velocity components toward C10 and C14, such as unresolved multiple sources, unresolved Keplerian disk or unresolved protostellar feedback within the condensation. Higher spatial and spectral resolution observations are required to distinguish these possibilities and determine the origin of these multiple velocity components.

Estimate of gas temperatures

We derived the gas temperatures (T_gas) using the K-ladder of CH₃CN J = 12–11 and ¹³CH₃CN J = 13–12 transitions with the XCLASS package⁵⁸. The Markov Chain Monte Carlo tasks built in XCLASS were used to explore the parameter space during the fitting process. For the combined data, the signal-to-noise ratios of the CH₃CN and ¹³CH₃CN are not sufficient to derive a reliable temperature map in the majority of cores, and they are not detected toward C12, C29 and three binary systems (C22–C38, C39–C42 and C35–C40). To improve the signal-to-noise ratio with minimal nearby source(s) contamination, we averaged the spectra within a half beam size toward the condensations. We exclude the K < 4 ladders in regions where these lower energy transitions become optically thick, that is, where line profile shows self-absorption or saturated emission. Although optical depth and chemical processes might affect the estimated temperatures, this is the best estimate of the temperature structure of the inner envelope we can obtain. To improve the fitting of the CH₃CN and ¹³CH₃CN lines, we include CH${}_{3}^{13}$CH J = 12–11 lines and other molecular lines (that is, CH₃OH, HNCO) in fitting for CH₃CN and CH₃OCH₃ for ¹³CH₃CN, if they are detected. The derived rotational temperatures range from 108 to 532 K (Table 1).

The rotational temperatures derived from ¹³CH₃CN are higher than those of CH₃CN. This is because the CH₃CN lines have a higher optical depth and preferably trace the surface of the structure, while ¹³CH₃CN is optically thinner and better trace the interior of the structure. This clearly suggests that these objects exhibit temperature gradients and are internally heated by the protostar(s) at the centre. Therefore, we use the temperature derived from ¹³CH₃CN to estimate the mass and luminosity, and in the case that ¹³CH₃CN is not sufficiently strong to allow a temperature measurement, we use the temperature derived from CH₃CN. Examples of spectra of ¹³CH₃CN and CH₃CN for the quadruple system are presented in Supplementary Fig. 1.

Computing the luminosity of the embedded protostar

With the derived gas temperature and taking into account the dust emissivity, we are able to approximately estimate the luminosity of the central heating source according to the relation between the temperature distribution and embedded protostellar luminosity^59,60, which is given by the following equation

$$L=1{0}^{5}\,{{\mbox{}}}{L}_{\odot }{{\mbox{}}}\,\times \,{\left(\frac{{T}_{\mathrm{D}}}{65{\,{\mbox{K}}}}\right)}^{4+\beta }{\left(\frac{0.1}{f}\right)}^{-1}{\left(\frac{0.1{\,{\mbox{pc}}}}{r}\right)}^{-2},$$

(2)

where T_D is the dust temperature at the radius r, β is the spectral index of dust emissivity at far-infrared wavelengths and f is its value at the wavelength of 50 μm. The β usually ranges from 0 to 1 and f is usually adopted as 0.1 (refs. ^60,61). The gas temperature derived from either ¹³CH₃CN or CH₃CN based on the averaged spectrum within a half of beam size of the condensation’s continuum peak can be used as a good approximation of T_D at the radius r = 130 AU (corresponding to the half beam size of ∼0.025″), where the densities are sufficiently high (>10^4.5 cm⁻³) for the dust and gas to be well coupled⁶². The ¹³CH₃CN and CH₃CN lines are not detected in C12, C29 and three binary systems (C22–C38, C39–C42, C35–C40). Thus, we refrain from estimating the M_* for these condensations to avoid the large uncertainty.

Assuming β = 1, the derived luminosities range from 50.9 to 1.5 × 10⁵ L_⊙ (Table 1), corresponding to spectral B9- to O6.5-type ZAMS stars³⁵, whose mass would be about 2.5–26.1 M_⊙ according to the mass–luminosity (M–L) relation³⁶. The total derived luminosity (∼4.2 × 10⁵ L_⊙) is higher than the value (∼2 × 10⁴ L_⊙) estimated on clump scale that has an uncertainty up to a factor of a few⁴⁵. We note that the derived temperatures could have some contaminations from neighbouring sources, and the luminosity of clump should be treated as a lower limit of total bolometric luminosity due to the lack of robust measurement at mid-infrared⁴⁵. The derived masses are regarded as conservative upper limits. There are five condensations (C1, C5, C8, C10 and C14) with estimated luminosities of 8.0 × 10³–1.5 × 10⁵ L_⊙, corresponding to a B1–O6.5 spectral type ZAMS star of >8 M_⊙ among the multiple systems (Table 1). The derived luminosity will be lower if a smaller β is adopted. However, even if β = 0 is adopted, there are still three condensations (C1, C10 and C14) associated with a M_* > 8 M_⊙. Therefore, a massive protostar should exist in both quintuple and quadruple systems, as also suggested by the presence of Class II CH₃OH maser, which are excited in high-density regions by strong radiation fields and exclusively tracing high-mass star-forming regions (Fig. 1).

Estimating ambient gas mass from dust continuum emission

The brightness temperatures of the dust emission in the condensations are lower than the gas temperatures T_gas, which is a good approximation of T_D. To check if the dust emission is optically thin, we computed the optical depth τ_cont of the continuum emission at the peak position of each condensation using⁶³

$${\tau }_{{{{\rm{\nu }}}}}^{{{{\rm{beam}}}}}=-{{{\rm{\ln }}}}\left(1\,-\,\frac{{S}_{{{{\rm{\nu }}}}}^{{{{\rm{beam}}}}}}{{{{\varOmega }}}_{{{{\rm{A}}}}}{B}_{\nu }({T}_{\mathrm{D}})}\right)$$

(3)

where B_ν is the Planck function at the dust temperature T_D, ${S}_{{{{\rm{\nu }}}}}^{{{{\rm{beam}}}}}$ is the continuum peak flux density and Ω_A is the beam solid angle. The condensations are dense enough (>10^4.5 cm⁻³) for gas and dust to be well coupled and in thermal equilibrium. As such, the gas temperature derived from the ¹³CH₃CN or CH₃CN should be approximately equal to the dust temperature. The derived optical depths are 0.04–0.27 with a mean value of 0.1 for all available condensations, indicating optically thin dust emission.

The observed 1.3 mm continuum emission is dominated by thermal dust emission toward G333 because the hydrogen recombination line (that is, H30α) is not detected toward condensations and the ATCA 3.3 mm continuum emission is also dominated by dust emission⁴¹. We calculate the ambient gas mass for the condensations following

$${M}_{{{{\rm{amb}}}}}=\eta \frac{{S}_{\nu }\,{{{{{D}}}^{2}}}}{{\kappa }_{\nu }\,{B}_{\nu }({T}_{\mathrm{D}})},$$

(4)

where η = 100 is the gas-to-dust ratio, S_ν is the measured integrated source flux, m_H is the mass of a hydrogen atom, μ = 2.8 is the mean molecular weight of the interstellar medium, D = 5.2 kpc is the distance to the source and κ_ν is the dust opacity at a frequency of ν. We adopted a value of 0.9 cm⁻² g⁻¹ for κ_1.3mm, which corresponds to the opacity of thin ice mantles and a gas density of 10⁶ cm⁻³ (ref. ⁶⁴). We use the lowest temperature of 108 K derived from CH₃CN as an approximation to the temperature for the condensations in which ¹³CH₃CN and CH₃CN are not detected. The actual temperature should be lower than 108 K, indicating that the derived mass is the lower limit for the condensations. The derived ambient gas masses of multiple systems are between 0.19 and 1.47 M_⊙ (Table 1), with the mean and median values of 0.59 and 0.52 M_⊙, respectively. The estimated ambient gas masses should be regarded as lower limits due to the interferometric observations suffering from missing flux.

Stability analysis of multiple system

To assess the stability of the multiple system, we compute the potential and kinetic energy of each object following the approach introduced in ref. ¹³. The gravitational potential energy, W_i, and kinetic energy, E_i, can be calculated by

$${W}_{i}=-\mathop{\sum}\limits_{i\ne j}\frac{G{m}_{i}{m}_{j}}{{r}_{ij}},$$

(5)

$${E}_{i}=\frac{1}{2}{m}_{i}{({{{{{\bf {V}}}}}}_{{{{\rm{i}}}}}-{{{{{\bf {V}}}}}}_{{{{\rm{com}}}}})}^{2},$$

(6)

where m_i and m_j are the masses of object i and j, r_ij is the separation between i and j, V_i is the (line-of-sight) velocity of object i, and V_com is the velocity of the centre of mass of the system. We determine the V_com through

$${{{{{\bf {V}}}}}}_{{{{\rm{com}}}}}=\frac{\mathop{\sum}\nolimits_{k}{m}_{k}{\mathbf{V}}_{k}}{\mathop{\sum}\nolimits_{k}{m}_{k}},$$

(7)

where m_k and V_k are the mass and velocity of the object k in the multiple system, respectively. A star with E_i/|W_i| < 1 is considered to be bound to the system.

The full velocity difference is $\sqrt{3}$ times the velocity difference along the line-of-sight, ΔV_3D = $\sqrt{3}\Delta {{{{{\bf{V}}}}}}_{{{{\rm{1D}}}}}$ = $\sqrt{3}$(V_i − V_com), assuming the measured velocity difference is representative of the one-dimensional velocity difference. Similarly, the total separation is $\sqrt{2}$ times the projected separation on the sky, r_3D = $\sqrt{2}{r}_{{{{\rm{1D}}}}}$ = $\sqrt{2}{r}_{ij}$, assuming that the measured projected separation is a good approximation of the separation along the line-of-sight. The observed mean separation, 〈r_1D〉, is about 731 AU, which is consistent with the typical projected value of 700 AU, r_1D = ${r}_{{{{\rm{3D}}}}}/\sqrt{2}$ = 1,000/$\sqrt{2}$ = 700 AU, in the simulation of multiple star formation via core fragmentation²⁹.

We used the ambient mass M_amb and protostar mass M_* to calculate the kinetic and gravitational energies for condensations in both one- and three-dimensional scenarios. We find that all multiple systems are gravitationally bound (Fig. 3), with exceptions for two condensations (C10 and C14) that have E_i/|W_i| > 1 for 3D velocity difference in the case of using M_amb. However, these two condensations are gravitationally bound if the central protostar mass is considered (Fig. 3). If the total mass, M_tot = M_am + M_*, is used, the E_i/|W_i| ratio will be smaller. Therefore, we conclude that all multiple systems are consistent with being gravitationally bound at the present stage. Table 1 presents the E_i and W_i for each condensation.