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

The dominant mechanism forming multiple stellar systems in the high-mass regime (M$_\ast \gtrsim $ 8 $M_{\odot}$) 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.

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.).The contour levels are [5, 8, 11, 14,17, 20, 23, 26, 29] × σ , where root mean squared (rms) noise of image is σ =0.6 mJy beam −1 .b, ALMA low-resolution (θ ∼ 0.32 ′′ ) 1.3 mm continuum image.The contour is the 7 σ , where σ =0.16 mJy beam −1 .The ellipses show the identified cores based on the low-resolution continuum image.The red pluses are the Class II CH 3 OH maser 29 .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 beam −1 .The dense cores (#1, #2, #3, ...) and condensations (C1, C2, C3, ...) are numbered in order of descending integrated intensity.c, Dense core #2 fragments into quadruple condensation system (C8, C10, C14, C17), and dense core #3 fragments into triple condensation system (C6, C12, C26).d, Dense core #1 fragments into quintuple condensation system (C1, C3, C4, C5, C16), and dense core #17 fragments into binary condensation (C11-C29).e, Dense core #7 fragments into binary condensation (C22-C38).f, Dense core #30 fragments into binary condensation (C39-C42).g, Dense core #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.
We report the direct imaging of 1 quintuple, 1 quadruple, 1 triple, and 4 binary systems in the high-mass protocluster G333.23-0.06(hereafter G333) by 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. 28 ) 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.
The projected separations of these multiple systems are between 327 and 1406 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 30 .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 ∼2100 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-ofsight 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 1664 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 3443 au (Table 1).The multiple condensation systems are embedded in the dense cores.Each condensation likely harbors 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, e.g., CH 3 CN, 13 CH 3 CN, and CH 3 OCHO (Fig. 2 and Extended Data Figs.5-7), and dense gas tracer, e.g., CH 3 OH, 13 CH 3 OH, SO 2 , SO, HC 3 N, HNCO, NH 2 CHO, H 2 CO, and H 13 2 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 3 au and significant different velocities 31 , is inconsistent with the observed separation distributions and small velocity differences.
If the parental disks are small and the separations of fragments are significantly widened within a short timescale 32 , we cannot completely rule out the possibility that disk fragmentation is responsible for the close binary systems (e.g., 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 5000-10000 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 33 , 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 instance, ongoing global collapse and dynamical interactions of multiple systems could lead to inward migration 34 , which moves the fragments closer together.

C14
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 ( √ 3 times the line-of-sight velocity difference) and on-sky separation (red symbols), (3) line-ofsight velocity difference and 3D separation ( √ 2 times the on-sky separation; blue symbols), (3) 3D velocity difference and 3D separation (orange symbols).The black dashed line marks E i /|W i | = 1.The members of quintuple system are marked with different color shadows, i.e., 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.

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 35 .
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 5 L ⊙ , corresponding to B9 to O6.5 spectral type ZAMS stars 36 .The masses are estimated to be between roughly 2.5 and 26.1 M ⊙ according to the mass-luminosity relation 37 (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 3 CN and its isotopologues.There are at least 3 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 3 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, etc.), 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-1400 au 38 , 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.

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.06G333.23-0.06 is a typical high-mass star-forming region [39][40][41][42] at a distance of 5.2 kpc 28 associated with Class II CH 3 OH maser emission 29 , 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 have a mass reservoir of ∼3000 M ⊙ with a mean column density of 1.6 × 10 23 cm −2 within a 1.2 pc radius 45 .As a representative example of high-mass star-forming region, G333 provides an ideal laboratory to study the binary and high-order multiple systems 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 shortbaseline or low-resolution) on November-05-2016 and 42 antennas in configuration similar to C43-8 (hereafter long-baseline) on July-28-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-1100 m and 91-8547 m for shortbaseline 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 −1 ) and a bandwidth of 1.875 GHz.These windows covered the frequency ranges of 233.5-235.5 GHz, 231.0-233.0GHz, 219.0-221.0GHz, 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 center used is (α(ICRS), δ (ICRS)) = 16h19m51.20s,−50 • 15 ′ 13. ′′ 00.
The visibility data calibration was performed using the CASA (version 5.4.0-70)software package 48 .We produced continuum data from line-free channels and continuum-subtracted data cubes for each observation epoch using the procedure described in ref. 46 .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 ′′ (1820 au×1560 au, panel b of Fig. 1) with a position angle of P.A. = -46.18• , 0.059 ′′ × 0.038 ′′ (307 au×198 au) with a P.A. = 56.23 • , and 0.066 ′′ × 0.039 ′′ (343 au×203 au, panels c-g of Fig. 1) with a P.A. = 54.47 • for short-baseline, long-baseline and combined dataset, respectively.The achieved 1σ rms noise levels continuum images are about 0.16 mJy beam −1 , 0.05 mJy beam −1 , and 0.05 mJy beam −1 for short-baseline, long-baseline, and combined data, respectively.Data cubes for each spectral window were produced using the automatic masking procedure YCLEAN 49 , which automatically cleans each map channel with custom-made masks.More details on the YCLEAN algorithm can be found in ref. 49 .The lines images 1σ rms noise are about 10 mJy beam −1 , 3 mJy beam −1 , and 3 mJy beam −1 with a channel width of ∼0.65 km s −1 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. 42 (panel a of Fig. 1).All images shown in this letter are prior to 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 3 − 10 4 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 50 .The performance of astrodendro in characterizing the dense structure parameters (e.g., 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 pre-select dense structures (i.e., the leaves in the terminology of astrodendro) from the 1.3 mm continuum images.We then use the parameters of the pre-select 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 axise 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 neighboring compact structures min − delta = 1σ ; 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 short-baseline 1.3 mm continuum image and 44 condensations in 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 3 OH 4 2,2 − 3 1,2 (E u /k = 45.46K) line that is detected in the majority of condensations in order 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), i.e., △V lsr < 1 km s −1 that is smaller than the line-of-sight velocity differences (2.0-9.5 km s −1 ) 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 3 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 51 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 6 -1 × 10 7 cm −3 , respectively.The derived thermal Jeans lengths range from 5000 to 10000 au.In the turbulent Jeans fragmentation scenario, the σ eff includes thermal and nonthermal velocity components, σ eff = c 2 s + σ 2 nth , under the assumption that the turbulence is acting as an isotropic of support.The measured line widths of CH 3 CN are 3-5 km s −1 for the parent cores, resulting turbulent Jeans lengths of 20000-54000 au.

Search for disk kinematic structure
The observations cover the typical disk tracers, including CH 3 CN and its isotopologues, and CH 3 OCHO, as well as other dense gas tracers, for instance H 2 CO and its isotopologues, CH 3 OH and its isotopologues, HC 3 N, NH 2 CHO, SO 2 , SO, HNCO, HCOOH, 13 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 which have a channel width of ∼0.65 km s −1 .The dense gas tracers are not sufficiently strong to allow a reliable determination of kinematic information for 3 binary systems (C22-C38, C39-C42, C35-C40) in the combined data.
There are 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 structure are found in the position-velocity (PV) diagram toward these cores.The velocity gradients trace either the outflows, or the large scale gas motions (e.g., gas flow, toroidal motions 52 ).These dense cores are associated with unipolar, bipolar, and/or perpendicular outflows identified by the SiO emission from the ALMA shortbaseline 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 53 .This further suggests that the quintuple, quadruple, and triple systems are formed from core fragmentation 54 .The detailed analysis of molecular outflows is beyond the scope of this letter, 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 (e.g., unresolved outflows or accretion flows).6][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 an 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 3 CN J = 12-11 and 13 CH 3 CN J = 13-12 transitions with the XCLASS package 58 .The Markov Chain Monte Carlo tasks built in XCLASS was used to explore the parameter space during the fitting process.For the com-bined data, the signal-to-noise ratios of the CH 3 CN and 13 CH 3 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, i.e., where line profile show 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 3 CN and 13 CH 3 CN lines, we include CH 13 3 CH J = 12-11 lines and other molecular lines (i.e., CH 3 OH, HNCO) in fitting for CH 3 CN and CH 3 OCH 3 for 13 CH 3 CN, if they are detected.The derived rotational temperatures range from 108 to 532 K (Table 1).
The rotational temperatures derived from 13 CH 3 CN are higher than those of CH 3 CN.This is because the CH 3 CN lines have a higher optical depth and preferably trace the surface of the structure, while 13 CH 3 CN is optically thinner and better trace the interior of the structure.This clearly suggests that these objects exhibit temperatures gradients and are internally heated by the protostar(s) at the centre.Therefore, we use the temperature derived from 13 CH 3 CN to estimate the mass and luminosity, and in the case that 13 CH 3 CN is not sufficiently strong to allow a temperature measurement, we use the temperature derived from CH 3 CN.Examples of spectra of 13 CH 3 CN and CH 3 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,61 , which is given by following equation 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 13 CH 3 CN or CH 3 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 −3 ) for the dust and gas to be well coupled 62 .The 13 CH 3 CN and CH 3 CN lines are not detected in C12, C29, and 3 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 5 L ⊙ (Table 1), corresponding to spectral B9-to O6.5-type ZAMS stars 36 , whose mass would be about 2.5 to 26.1 M ⊙ according to the mass-luminosity (M-L) relation 37 .The total derived luminosity (∼4.2×10 5 L ⊙ ) is higher than the value (∼2×10 4 L ⊙ ) estimated on clump scale that has an uncertainty up to a factor of a few 45 .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 midinfrared 45 .The derived masses are regarded as conservative upper limits.There are 5 condensations (C1, C5, C8, C10, and C14) with estimated luminosities of 8.0×10 3 -1.5×10 5 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 3 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 3 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 63 where B ν is the Planck function at the dust temperature T D , S beam ν is the continuum peak flux density, Ω A is the beam solid angle.The condensations are dense enough (>10 4.5 cm −3 ) for gas and dust to be well coupled and in thermal equilibrium.As such, the gas temperature derived from the 13 CH 3 CN or CH 3 CN should be approximately equal to the dust temperature.The derived optical depths are 0.04-0.27with 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 (i.e., H30α) is not detected toward condensations and the ATCA 3.3 mm continuum emission is also dominated by dust emission 42 .We calculate the ambient gas mass for the condensations following where η = 100 is the gas-to-dust ratio, S ν is the measured integrated source flux, m H is the mass of an 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 −2 g −1 for κ 1.3mm , which corresponds to the opacity of thin ice mantles and a gas density of 10 6 cm −3 (ref. 64).We use the lowest temperature of 108 K derived from CH 3 CN as an approximation to the temperature for the condensations in which 13 CH 3 CN and CH 3 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. 13 .The gravitational potential energy, W i , and kinetic energy, E i , can be caculated by where m i and m j are the masses of object i and j, r i j 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 where m k and V k are the mass and velocity of the object k in the multiple system.A star with E i /|W i | < 1 is considered to be bound to the system.The full velocity difference is √ 3 times the velocity difference along the line-of-sight, , assuming the measured velocity difference is representative of the one-dimensional velocity difference.Similarly, the total separation is √ 2 times the projected separation on the sky, r 3D = √ 2r 1D = √ 2r i j , 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 3D / √ 2 = 1000/ √ 2 = 700 au, in the simulation of multiple star formation via core fragmentation 30 .
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 amb + 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.

Fig. 2 |
Fig. 2 | Examples of intensity-weighted velocity maps and position-velocity (PV) diagrams of molecular lines for the quintuple system.Columns left, middle, and right present the results derived from the CH 3 CN 12 4 − 11 4 , 13 CH 3 CN 13 3 − 12 3 , and CH 3 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 3 CN 12 4 − 11 4 (a), 13 CH 3 CN 13 3 − 12 3 (b), and CH 3 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 3 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 is 7.5 mJy beam −1 km s −1 .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 3 CN 12 4 − 11 4 (d), 13 CH 3 CN 13 3 − 12 3 (e), and CH 3 OCHO 20 0,20 − 19 0,19 (f).Contours levels start at 4σ rms and increase in step of 1σ rms interval, where σ rms is 1.2 mJy beam −1 .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 the dense core #1 and its systemic velocity of V lsr = -85 km s −1 .

. 1 |. 3 |. 5 |. 6 |b 8 |
ALMA high-resolution 1.3 mm continuum image.The cyan ellipses are the identified dense cores as shown in Extended Data Fig. 4. The green crosses show condensations identified from ALMA high-resolution 1.3 mm continuum image.The grey contour shows the 7σ , where σ = 0.05 mJy beam −1 .The synthesized beam size of 1.3 mm continuum image present in the lower left corner with a white ellipse.Extended Data Fig. 2 | The separation distributions of multiple systems.The separations range from 327 to 1406 au, with a mean value of 731 au.Intensity-weighted velocity maps of CH 3 OH derived from the high-resolution data toward multiple systems.The grey contours show the corresponding velocity-integrated intensity maps.Contours levels start at 3σ rms and increase in step of 3σ rms interval, where σ rms is 7.6, 6, 6.4, 6.4, 6, and 5 mJy beam −1 km s −1 for a-f.The condensation C40 does not have sufficient signal-to-noise (SNR) ratio for the moment maps, while the condensation averaged spectrum has sufficient SNR to determine the centroid velocity.The black ellipses in the lower left corner of each panel denote the synthesized beam of lines images.Intensity-weighted velocity maps derived from the low-resolution data toward parent cores of multiple systems.a-f, We show intensity-weighted velocity maps of CH 3 CN 12 4 − 11 4 (a and d), 13 CH 3 CN 13 3 − 12 3 (b and e), and CH 3 OCHO 20 0,20 − 19 0,19 (c and f) for dense cores #1, #17, #2, and #3.The blue and red arrows show the directions of the outflows seen in the SiO emission from the ALMA low-resolution data.g-i, Intensity-weighted velocity maps of H 2 CO 3 2,2 − 2 2,1 for dense cores #7 (g) and #15 (i).The magenta ellipses and black crosses show the dense cores and their embedded condensations, respectively.The red plusses marks the Class II CH 3 OH maser positions.The grey contours show the velocity-integrated intensity maps with levels of [5, 10, 15, 20, 40, 60, 80, 100]×σ rms , where the σ rms is 13 mJy beam −1 km s −1 for a-f and 9 mJy beam −1 km s −1 for g-i.The black ellipses in the lower left corner of each panel denote the synthesized beam of lines images.Intensity-weighted velocity maps derived from the high-resolution data toward multiple systems.a-c, Intensity-weighted velocity maps of CH 3 CN 12 4 − 11 4 (a), 13 CH 3 CN 13 3 − 12 3 (b), and CH 3 OCHO 20 0,20 − 19 0,19 (c) for the quadruple system.d-e, Intensity-weighted velocity maps of CH 3 OH 4 2,2 − 3 1,2 (d) and SO 5, 6 − 4, 5 (e) for the triple and the binary systems, respectively.The black and magenta ellipses show the condensations and their parent cores, respectively.The red plusses marks the Class II CH 3 OH maser positions.The grey contours show the velocity-integrated intensity maps with levels of [5, 10, 15, 20, 30]×σ rms , where the σ rms is 9 mJy beam −1 km s −1 for a-c and 6 mJy beam −1 km s −1 for d-e.The black ellipses in the lower left corner of each panel denote the synthesized beam of lines images.a Extended Data Fig. 7 | Intensity-weighted velocity maps of CH 3 CN 12 4 − 11 4 for two velocity ranges of [-98, -91] km s −1 and [-90, -84] km s −1 .The black and magenta ellipses show the condensations and their parent cores, respectively.The red plusses marks the Class II CH 3 OH maser positions.The grey contours show the velocity-integrated intensity maps with levels of [5, 10, 15, 20, 30]×σ rms , where the σ rms is 6 mJy beam −1 km s −1 .The black ellipses in the lower left corner of each panel denote the synthesized beam of lines images.Extended Data Fig.Molecular outflows seen in SiO emission in the ALMA lowresolution data.SiO redshifted and blueshifted contours overlaid on the ALMA low-resolution 1.3 mm continuum.The green and magenta arrows present the blueshifted and redshifted directions of the SiO outflow.The cyan ellipses show the dense cores.In panel a, the SiO emission shows a prominent blueshifted, a clear redshifted, and a weak redshifted component from dense core #1.In panel b, the SiO emission reveals complicated outflow spatial morphologies toward dense cores #2 and #3.Dense core #3 drives two blueshifted outflows.Dense core #2 drives at least two blueshifted and two redshifted outflows.Contours levels start at 3 σ rms and increase in step of 1.5 σ rms interval, where σ rms is 0.023 Jy beam −1 km s −1 .The synthesized beam size of 1.3 mm continuum image present in the lower left corner.The detected misaligned outflows further suggest that the quintuple, quadruple, and triple systems are formed from core fragmentation 54 .

Table 1 |
Properties of condensations