The ALMA-QUARKS survey: A Keplerian Circumbinary Disk Candidate in the High-mass Star-forming Region IRAS 08303-4303

Email: [email protected], [email protected]

Introduction

High-mass stars ($M \ge8\ M_\odot$) typically form within dense regions of molecular clouds, known as hubs, where filaments converge (S. Dib et al. 2020; M. S. N. Kumar et al. 2020; J. W. Zhou et al. 2025). These potential high-mass star-forming regions may exist as prestellar cores for $< 10^4$ years (A. Duarte-Cabral et al. 2013; J. Tigé et al. 2017; F. Motte et al. 2018). During the early stages of high-mass star formation, stellar embryos with masses less than $8\ M_\odot$ are IR-quiet protostars (H. Beuther et al. 2007; F. Motte et al. 2007, 2018; D. Russeil et al. 2010). Protostars become IR-bright when the stellar embryos reach $M \ge8\ M_\odot$, entering the high-mass protostellar object (HMPO) stage (F. Motte et al. 2007, 2018; D. Russeil et al. 2010). The HMPO stage can be further divided based on chemical characteristics into two categories: chemically poorer (earlier) HMPOs (T. Gerner et al. 2014) and later, hot molecular cores (HMCs) that exhibit molecular-rich chemistry (T. Gerner et al. 2014; H.-L. Liu et al. 2021; S.-L. Qin et al. 2022). HMCs are characterized by rich spectra of complex organic molecules (COMs), high temperatures ($\gtrsim100$ K), and compact source sizes ($\lesssim0.1$ pc) (R. Cesaroni 2005). Eventually, when the ultraviolet (UV) radiation from the embedded high-mass (proto)star ionizes the surrounding gas, a hyper-compact HII (HCHII) or ultra-compact HII (UCHII) region is formed. Due to the short evolutionary timescales, multiple stages can overlap; for instance, line-rich HMC chemistry and emission from ionized gas may coexist within the same core (H.-L. Liu et al. 2021). Additionally, because high-mass stars typically form in clusters, individual cores at varying evolutionary stages are frequently observed in close proximity to one another (T. Liu et al. 2020).

Around an accreting protostar, the centrifugal force of rotating material typically balances gravity, leading to the formation of a rotationally supported structure. While this often results in a Keplerian disk ($v \propto r^{-0.5}$, B. Zhao et al. 2020), actual disks frequently exhibit deviations from a strict Keplerian profile due to additional physical factors, such as self-gravity, internal gas pressure in hot disks, and magnetic braking effects. Some of the infalling material is ejected as a bipolar outflow, a process driven by complex magnetohydrodynamics (MHD) that helps carry away excess angular momentum (K. Tomisaka 2000). This disk-outflow system is crucial for massive young stellar objects (MYSOs), as it protects the accretion process from the intense radiation ($L \ge10^3\ L_\odot$) pressure exerted by the central protostars (R. Kuiper et al. 2010). Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) has identified Keplerian disk candidates around several MYSOs (M. T. Beltrán et al. 2014; K. G. Johnston et al. 2015; J. D. Ilee et al. 2016, 2018; A. Ginsburg et al. 2018; L. T. Maud et al. 2019; L. Moscadelli et al. 2019; X. Lu et al. 2022; G. M. Williams et al. 2022), some of which extend up to a few hundred au. These disks are often associated with COMs such as acetonitrile (CH$_3$CN, e.g., K. G. Johnston et al. 2015; J. D. Ilee et al. 2016) and methanol (CH$_3$OH, e.g., J. D. Ilee et al. 2016; G. M. Williams et al. 2022). In certain cases, refractory molecules like sodium chloride (NaCl, e.g., A. Ginsburg et al. 2019, 2023; K. E. I. Tanaka et al. 2020) or hydrogen recombination lines from ionized gas (e.g., R. Miyawaki et al. 2023) can be used to trace the presence of these disks. In contrast, bipolar molecular outflows are usually traced using carbon monoxide (CO, e.g., A. Y. Yang et al. 2022). Silicon monoxide (SiO) line emissions frequently indicate central jet-like features along the axes of broader CO outflows (e.g., J. A. López-Vázquez et al. 2024).

High-mass stars are often found in binary or multiple star systems. Multiple surveys, mainly in the optical and near-infrared regimes, have shown that the multiplicity fraction—defined as the proportion of primaries having at least one companion—among main-sequence high-mass stars exceeds 90% (S. S. R. Offner et al. 2023). Despite this high percentage, the early formation phases of such massive binary systems remain largely unexplored. In low-mass binary systems, we can observe large rotating structures that surround both protostars and their accretion disks, commonly referred to as “circumbinary disks”, e.g., L1551 NE (S. Takakuwa et al. 2012, 2014), L1551 IRS 5 (T.-L. Chou et al. 2014), [BHB 2007] 11 (F. O. Alves et al. 2019). However, circumbinary disks have only been observed around a very limited number of massive proto-binary systems. One example is the massive protobinary system IRAS 16547–4247 in which a circumbinary disk was detected using the well-known disk tracer $\mathrm{CH_3CN}$ (K. E. I. Tanaka et al. 2020), along with other COMs, such as methyl formate ($\mathrm{CH_3OCHO}$, L. A. Zapata et al. 2019).

MYSOs are rare and are predominantly found at distances exceeding 2 kpc. This distance poses significant challenges for observational studies due to limitations in spatial resolution. However, the deployment of the ALMA has made it possible to resolve the small structures surrounding these distant objects. One notable statistical survey focused on high-mass star-forming regions in the galactic plane is the “ALMA Three-millimeter Observations of Massive Star-forming regions (ATOMS)” project (T. Liu et al. 2020). The ATOMS project observed 146 massive protoclusters at a wavelength of 3 mm with resolutions ranging from 1″ to 2″. Additionally, 139 of these protoclusters were further observed at a wavelength of 1.3 mm by the “Querying Underlying mechanisms of massive star formation with ALMA-Resolved gas Kinematics and Structures (QUARKS)” project (X. Liu et al. 2024). The QUARKS project achieved high resolutions of $\sim$0.″3, while also covering large-scale flux up to 27″ through Atacama Compact 7-m Array (ACA) observations (F. Xu et al. 2024). The QUARKS frequency band captures numerous molecular tracers of disks and outflows (as detailed in Table 2 of X. Liu et al. 2024), making it particularly well-suited for investigating various aspects of disk-outflow systems surrounding MYSOs.

Line	Frequency$^{a}$ (GHz)	$E_u^{a}$ (K)	NE	C1	C2	SW1$^{e}$	SW2	E$^{c}$
$\mathrm{CH_3CN}\ J_K = 12_0$–$11_0$	220.7472617	68.9	Y	Y	Y	Y	Y	Y$^{c}$
$\mathrm{CH_3CN}\ J_K = 12_1$–$11_1$	220.7430111	76.0	n	Y	Y	Y$^{c}$	Y	Y
$\mathrm{CH_3CN}\ J_K = 12_2$–$11_2$	220.7302611	97.4	n	Y	Y	n	Y	n
$\mathrm{CH_3CN}\ J_K = 12_3$–$11_3$	220.7090170	133.2	Y$^{c}$	Y	Y	Y$^{c}$	Y	n
$\mathrm{CH_3CN}\ J_K = 12_4$–$11_4$	220.6792874	183.1	n	Y	Y	n	Y	n
$\mathrm{CH_3CN}\ J_K = 12_5$–$11_5$	220.6410844	247.4	n	Y	Y	n	Y	n
$\mathrm{CH_3CN}\ J_K = 12_6$–$11_6$	220.5944237	325.9	n	Y	Y	n	Y	n
$\mathrm{CH_3CN}\ J_K = 12_7$–$11_7$	220.5393240	418.6	n	Y	Y	n	n	n
$\mathrm{CH_3CN}\ J_K = 12_8$–$11_8$	220.4758078	525.6	n	Y	Y$^{c}$	n	n	n
$\mathrm{CO}\ J = 2$–$1$	230.5380000	16.6	Y	Y	Y	Y	Y	Y
${}^{13}\mathrm{CO}\ J = 2$–$1$	220.3986842	15.9	Y	Y	Y	Y	Y	Y
$\mathrm{C}^{18}\mathrm{O}\ J = 2$–$1$	219.5603577	15.8	Y	Y	Y	Y	Y	Y
$\mathrm{SiO}\ J = 5$–$4$	217.104919	31.3	Y$^{d}$	Y	Y	n	Y	n
$\mathrm{H_{30}\alpha}$	231.900928$^{b}$	–	n	n	n	n	n	n
$\mathrm{SO}\ J_N = 6_5$–$5_4$	219.94944	35.0	Y	Y	Y	Y	Y	Y
$\mathrm{SO_2}\ J_{K_a,K_c} = 28_{3,25}$–$28_{2,26}$	234.187057	403.0	n	Y	Y	n	Y$^{c}$	n
$\mathrm{SO_2}\ J_{K_a,K_c} = 16_{6,10}$–$17_{5,13}$	234.4215884	213.3	n	Y	Y	n	Y$^{c}$	n
$\mathrm{HC_3N}\ J = 24$–$23$	218.32472	131.0	Y$^{c}$	Y	Y	Y	Y	Y$^{c}$
$\mathrm{H_2CO}\ J_{K_a,K_c} = 3_{0,3}$–$2_{0,2}$	218.22219	21.0	Y	Y	Y	Y	Y	Y
$\mathrm{H_2CO}\ J_{K_a,K_c} = 3_{2,2}$–$2_{2,1}$	218.47563	68.1	Y	Y	Y	Y	Y	Y
$\mathrm{H_2CO}\ J_{K_a,K_c} = 3_{2,1}$–$2_{2,0}$	218.76007	68.1	Y	Y	Y	Y	Y	Y
$\mathrm{N_2D^+}\ J = 3$–$2$	231.321828	22.2	n	n	n	n	n	n
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 4_{2,3}$–$3_{1,2}\ E$	218.44006	45.5	Y	Y	Y	Y	Y	Y
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 8_{0,8}$–$7_{1,6}\ E$	220.078561	96.6	Y$^{c}$	Y	Y	Y	Y	Y
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 10_{2,9}$–$9_{3,6}\ A^-$	231.28111	165.3	Y$^{c}$	Y	Y	Y$^{c}$	Y	Y$^{c}$
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 10_{-3,7}$–$11_{-2,9}\ E$	232.94580	190.4	n	Y	Y	n	Y	n
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 6_{1,5}$–$7_{2,5}\ A^-$, $v_t = 1$	217.29921	373.9	Y$^{c}$	Y	Y	n	Y	n
$\mathrm{CH_3OH}\ J_{K_a,K_c} = 18_{3,16}$–$17_{4,13}\ A^+$	232.78345	446.5	n	Y	Y	n	Y	n
$\mathrm{DCN}\ J = 3$–$2$	217.2385378	20.9	Y	Y	Y	Y	Y	Y
$\mathrm{H^{13}CO}\ J_{K_a,K_c} = 3_{1,2}$–$2_{1,1}$	219.90852	32.9	n	Y	Y	n	Y	Y
$\mathrm{HNCO}\ J_{K_a,K_c} = 10_{0,10}$–$9_{0,9}$	219.798274	58.0	n	Y	Y	n	Y	Y$^{c}$
$\mathrm{CH_3OCH_3}\ J_{K_a,K_c} = 13_{0,13}$–$12_{1,12}$	231.987858$^{f}$	80.9	n	Y	Y	Y$^{c}$	Y	Y$^{c}$
${}^{13}\mathrm{CS}\ J = 5$–$4$	231.220685	33.3	Y	Y	Y	Y	Y	n
$\mathrm{OCS}\ J = 19$–$18$	231.060993	110.9	Y	Y	Y	Y	Y	n

IRAS 08303–4303 is one of the targets in the ATOMS and QUARKS projects, located at a distance of 2.35 ± 0.55 kpc (H.-L. Liu et al. 2021; X. Liu et al. 2024). S. Faúndez et al. (2004) estimated the overall clump mass of this object to be 260 $\mathrm{M_\odot}$, based on 1.2 mm continuum flux measurements obtained with the Swedish-ESO submillimeter telescope. They also derived the bolometric luminosity to be $\sim$ 7000 $\mathrm{L_\odot}$ from the spectral energy distribution (SED) fitting, which incorporated both the 1.2 mm continuum flux and the fluxes from the four IRAS bands. X. Liu et al. (2024) determined the systemic LSR (Local Standard of Rest) velocity ($v_{\mathrm{sys}}$) of the overall clump to be approximately 14.5 km s$^{-1}$, based on $\mathrm{H^{13}CO^+}$ ATOMS data. The ATOMS project identified more than 20 emission lines towards the brightest core. Many of these emissions are from COMs (Z.-Y. Li et al. 2025a). The calculated excitation temperature using $\mathrm{CH_3OH}$, $\mathrm{CH_3OCHO}$, and propionitrile ($\mathrm{C_2H_5CN}$) is $\sim$100 K, which supports the hot core chemistry (S.-L. Qin et al. 2022). It is worth noting that the $\mathrm{H_{40\alpha}}$ hydrogen recombination line emission, which traces ionized HII regions, was not detected. Consequently, H.-L. Liu et al. (2021) classified it as a pure and strong hot molecular core candidate (pure s-cHMC). This classification is based on its high line richness ($N_{\mathrm{line}} \geq20$) and the non-detection of the $\mathrm{H_{40\alpha}}$ line. However, the ATOMS observations, with a resolution of $\sim$1.″5 (about 3500 au at a distance of 2.35 kpc), have not resolved the IRAS 08303–4303 region down to the scale of the disk.

In this study, we investigate the physical and dynamical properties of dense condensations and disk-outflow systems in the IRAS 08303–4303 region by primarily analyzing the QUARKS data at a resolution of $\sim$0.″3 (about 700 au). This level of detail has enabled us to identify a candidate for a circumbinary disk within the aforementioned hot core. While recent higher-resolution ($\sim$ 0.″05) observations have further resolved the dense condensations down to the $\sim$100 au scale (K. Yang et al. 2026), our analysis provides a crucial link between the core-scale (ATOMS) and these small-scale structures. The observational details are summarized in § 2, while the analysis results are presented in § 3. § 4 examines the kinematics of the circumbinary disk and its associated distinct features, and § 5 summarizes the main findings of our research.

Observation

ALMA observations of the massive protostellar system IRAS 08303–4303 were conducted in Bands 3 and 6. The Band 3 observations took place on 2019 October 14 and 19, as part of the ATOMS program (Project ID: 2019.1.00685.S; PI: Tie Liu). The first observation was conducted with the ACA, and the second with the 12-m array in the C-4 configuration. Among the eight spectral windows (SPWs) used, the continuum image was generated from two wide SPWs, each with a bandwidth of 1875 MHz. The combined Band 3 continuum image from the ACA and 12-m array has a synthesized beam of 1\farcs55 \times1\farcs16 with a position angle (PA) of -71\fdg3 and a noise level of 90 $\mu$Jy beam$^{-1}$. Detailed information regarding the ATOMS observations can be found in T. Liu et al. (2020).

The Band 6 observations were carried out on 2021 November 6, 2022 December 29, and 2023 April 8, as part of the QUARKS program (Project IDs: 2021.1.00095.S, 2022.1.00298.S, and 2023.1.00425.S; PI: Lei Zhu). The observations utilized three configurations: the ACA and the 12-m array in C-3 and C-5 configurations. A comprehensive overview of each configuration’s observation details is presented in Table 1. The three configuration data sets were calibrated using the pipeline of the Common Astronomy Software Applications (CASA; CASA Team et al. (2022)) and subsequently combined with the CASA task concat before imaging. For continuum imaging, we utilized all four 1875 MHz wideband spectral windows (SPWs) centered at 217.9, 220.3, 231.3, and 233.5 GHz, after manually masking the spectral lines. For both line and continuum imaging, we set the robust parameter to +0.5 in the CASA tclean task. Deconvolution was performed using the ‘multiscale’ algorithm for line imaging and the ‘mtmfs’ (Multi-Term Multi-Frequency Synthesis) algorithm for continuum imaging. Following this, we performed self-calibration in both phase and amplitude using the resulting 1.3 mm continuum image as a model, applying a solution interval equal to the length of a scan ($\sim60$ sec). We then applied these self-calibration solutions to the spectral data. Finally, we generated line cubes for all four SPWs after subtracting the continuum using the CASA task uvcontsub. The combined Band 6 continuum image resulted in a synthesized beam of 0\farcs33 \times0\farcs33 (PA = +66\fdg1) and a noise level of 0.15 mJy beam$^{-1}$. All spectral line images have synthesized beams between 0\farcs3 and 0\farcs4, with a velocity resolution of $\sim1.3$ km s$^{-1}$ and a noise level of 3 mJy beam$^{-1}$. These four SPWs were specifically designed to target various spectral lines, including CO, SiO, CH$_3$OH, and CH$_3$CN. Detailed properties and detections of these targeted lines are presented in Table 2. The maximum recoverable scale for this ALMA observation is 20″ ($\sim0.2$ pc). More details about the QUARKS observations can be found in X. Liu et al. (2024).

Configuration	Date	BL (m)	AR ($^{\prime\prime}$)	MRS ($^{\prime\prime}$)	N Ant.	ToS (min)	Calibrators
ACA	2021-11-06	8.9–45.0	4.9	20.0	9	5	J1058+0133, J0922–3959
C-3	2022-12-29	15.1–500.6	0.7	7.3	46	1	J1037–2934, J0922–3959
C-5	2023-04-08	15.3–1397.4	0.2	3.4	44	4	J1107–4449, J0922–3959, J0849–3541

Results

Dense Continuum Condensations

Figure 1 presents the continuum emission images of ATOMS and QUARKS. A central hot core in the lower-resolution 3 mm image ((S.-L. Qin et al. 2022)) is resolved into two distinct condensations, labeled ‘C1’ and ‘C2,’ when viewed in the higher-resolution 1.3 mm image. The projected distance between C1 and C2 is only 0\farcs38 ($\sim900$ au), indicating that they likely form a binary system. Notably, C2 emits stronger continuum emission than C1 and is considered the primary component. Six condensations have been identified visually; five bright and dense condensations labeled NE, C1, C2, SW1, and SW2, and a faint and isolated condensation labeled E (Figure 1(b)). All six condensations are aligned along a filamentary structure elongated in the northeast-southwest direction.

Figure 2 presents the spectra extracted at the continuum peaks of the six dense condensations. Both C1 and C2 display line-rich HMC spectra, including multiple transitions of CH$_3$CN and CH$_3$OH. SW2 also shows these transitions, but their intensities are much weaker. In contrast, NE, SW1, and E show line-poor spectra, displaying only a few detectable lines, such as transitions of CO isotopologues and H$_2$CO, with only marginal detections otherwise. These transitions exhibit extended morphologies, which means that some may be affected by contamination from the nearby hot cores C1 and C2. This suggests that these condensations might be at a colder, earlier evolutionary stage than C1, C2, and SW2. Alternatively, they could be lower-mass objects. The NE condensation contains a stellar embryo at its center, as strong bipolar outflows have been detected (see Figure 3). Moreover, its outflow parameters are significantly larger than those of typical low-mass YSOs and comparable to those of MYSOs (see § 3.2). Thus, this line-poor signature implies that NE could be classified as a line-poor HMPO prior to the formation of a hot core, rather than as a low-mass YSO. However, there is no clear sign of star formation activity in the SW1 and E condensations.

The H$_{30\alpha}$ line emission was not detected in any of the six condensations, indicating that none of them have reached the UCHII phase. As mentioned above, the H$_{40\alpha}$ line emission was not detected, either ((H.-L. Liu et al. 2021)). While an emission line was detected in C1 near the rest frequency of the H$_{30\alpha}$, its linewidth (full width at half maximum, FWHM), 3.5 km s$^{-1}$, differs greatly from the broad linewidth of 25–30 km s$^{-1}$ typically expected in UCHII regions ((D. O. S. Wood & E. Churchwell 1989; W. J. Kim et al. 2017, 2018)). Thus, it is likely one of the $^{33}$SO$_2$, CH$_3$OCH$_3$, or C$_2$H$_5$CN lines, which have rest frequencies within 1 MHz of H30$\alpha$. However, since there is no cm-wavelength data available, the presence of an HC/UCH$_{\rm II}$ region cannot be completely ruled out. Table 2 summarizes representative molecular lines and their detection toward each condensation. These lines were identified by comparing observed frequencies with those listed in the Cologne Database for Molecular Spectroscopy (CDMS, H. S. P. Müller et al. 2001, 2005; C. P. Endres et al. 2016).

The source position, deconvolved size (FWHM), position angle (PA), peak intensity, and integrated flux density of six dense condensations were derived from two-dimensional Gaussian fits using the CASA task imfit.

Five additional diffuse, extended components were also fitted along with the dense condensations to perform proper fittings around bright peaks after subtracting the surrounding diffuse emission. The temperatures of the C1, C2, SW2 condensations and the central diffuse component were determined from the rotation diagram analyses of multiple CH$_3$OH transitions, assuming the case of optically thin line emission in the local thermal equilibrium (LTE) condition (see Appendix § 7 for details). The dust temperature ($T_{\rm d}$) was assumed to be equal to the derived rotational temperature. For the other condensations, the CH$_3$OH transitions with $E_u > 100$ K were not significantly detected, so the temperatures cannot be calculated from rotational dia-

$^{a}$Deconvolved quantities

$^{b}$Temperature was calculated using CH$_3$OH (see Appendix § 7)

gram analysis. Three formaldehyde (H$_2$CO) transitions were detected in all dense condensations, but the H$_2$CO $3_{0,3}$–$2_{0,2}$ line is too optically thick to apply the method of X. D. Tang et al. (2017, 2018).

For the NE condensation, the line-poor spectrum and strong bipolar outflow suggest that it is a line-poor HMPO. The typical temperature for such an object is 50 K (T. Gerner et al. 2014), but the continuum condensation observed in this study has a much smaller scale than the ones detected in previous single-dish observations, so the temperature could be higher. However, it is expected to be lower than the typical hot core temperature of $\gtrsim100$ K, giving a range of 50–100 K. Assuming a uniform distribution of values between 50 and 100 K, this can be represented as $75 \pm14$ K (the standard deviation of a uniform distribution on $[a,b]$ is $(b-a)/\sqrt{12}$). For the SW1 and E condensations, star formation activity is unclear, and the temperature is likely lower than the 50 K typical of HMPOs. However, it might be higher than about 25 K because N$_2$D$^+$ was not detected toward either condensation. N$_2$D$^+$ is known to be depleted above the CO sublimation temperature ($\sim25$ K) due to reactions with CO and electrons (S. Li et al. 2022; T. Sakai et al. 2022). Similarly to NE, assuming a uniform distribution of temperatures between 25 and 50 K, this can be represented as $38 \pm7$ K.

The masses of the condensations were calculated using the following equation (H.-L. Liu et al. 2021) assuming that the 1.3 mm continuum is optically thin dust emission,

Here, $d$ is the distance to the source, $F_\nu$ is the integrated flux density (at $\nu= 225.7$ GHz) over the source size, $R_{\mathrm{gd}}$ is the gas-to-dust mass ratio (assumed to be 100), $B_\nu(T_{\mathrm{d}})$ is the Planck function at temperature $T_{\mathrm{d}}$, and $\kappa_\nu$ is the dust absorption coefficient, assumed to be 0.9 cm$^2$ g$^{-1}$ at 1.3 mm. This value corresponds to grains covered by a thin ice mantle and coagulated at gas densities of $10^6$ cm$^{-3}$ over a coagulation timescale of $10^5$ years (V. Ossenkopf & T. Henning 1994). We adopt uncertainties of 23% for $R_{\mathrm{gd}}$ and 28% for $\kappa_\nu$ in the calculation (P. S. Shanchez et al. 2017, 2019).

All the dense condensations have masses of $M \lesssim1 M_\odot$ (Table 3). It should be noted that these are not the protostellar masses but the masses of the surrounding dense gases; much of the mass may have already been accreted onto the central objects. Additionally, a significant fraction of the cold, lower-density gas may not be observed owing to sensitivity limits or being resolved out, and therefore not included in the mass calculation. There is also a diffuse continuum component surrounding the central binary system. Using the same assumptions as before, the temperature of this component is estimated to be $127 \pm2$ K based on CH$_3$OH, which corresponds to a component mass of $\sim2 M_\odot$. This component may represent the circumbinary disk around the protobinary system (see § 3.3 and § 4.1 for further details). The observed properties of the six dense condensations and the surrounding diffuse components are summarized in Table 3.

Name	Center (J2000)		FWHM$^{a}$	PA$^{a}$	$I_{\mathrm{peak}}$	$F_{1.3\ \mathrm{mm}}$	Temperature	Mass
	R.A.	Dec.	($'' \times''$)	($^\circ$)	(mJy beam$^{-1}$)	(mJy)	(K)	($M_\odot$)
NE	08$^\mathrm{h}$32$^\mathrm{m}$09$^\mathrm{s}$.05	$-43^\circ13'43.''4$	$0.50 \pm0.02 \times0.32 \pm0.02$	$95 \pm4$	$12.6 \pm0.3$	$31.8 \pm1.0$	$75 \pm14$	$0.9 \pm0.5$
C1	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.66	$-43^\circ13'45.''4$	$0.30 \pm0.01 \times0.22 \pm0.01$	$60 \pm7$	$14.9 \pm0.3$	$24.4 \pm0.7$	$235 \pm3^{b}$	$0.2 \pm0.1$
C2	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.64	$-43^\circ13'45.''9$	$0.55 \pm0.01 \times0.22 \pm0.01$	$147 \pm1$	$19.6 \pm0.3$	$46.0 \pm0.9$	$171 \pm6^{b}$	$0.5 \pm0.3$
SW1	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.46	$-43^\circ13'47.''3$	$0.11 \pm0.04 \times0.08 \pm0.04$	$166 \pm115$	$13.3 \pm0.3$	$14.4 \pm0.5$	$38 \pm7$	$0.8 \pm0.5$
SW2	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.45	$-43^\circ13'49.''3$	$0.69 \pm0.05 \times0.21 \pm0.03$	$67 \pm2$	$6.4 \pm0.3$	$17.2 \pm1.0$	$262 \pm9^{b}$	$0.13 \pm0.06$
E	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.79	$-43^\circ13'46.''3$	$0.28 \pm0.07 \times0.22 \pm0.08$	$57 \pm40$	$3.2 \pm0.3$	$5.1 \pm0.7$	$38 \pm7$	$0.3 \pm0.2$
Diffuse Components
	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.97	$-43^\circ13'41.''9$	$6.2 \pm0.3 \times2.3 \pm0.1$	$80 \pm2$	$1.29 \pm0.07$	$170 \pm9$
Central	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.65	$-43^\circ13'45.''5$	$2.59 \pm0.09 \times1.59 \pm0.06$	$16 \pm3$	$3.5 \pm0.1$	$136 \pm5$	$127 \pm2^{b}$	$2.1 \pm1.0$
	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.42	$-43^\circ13'48.''5$	$4.9 \pm0.3 \times1.5 \pm0.1$	$31 \pm2$	$1.5 \pm0.1$	$106 \pm6$
	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.49	$-43^\circ13'48.''8$	$1.80 \pm0.10 \times0.36 \pm0.03$	$14 \pm1$	$4.5 \pm0.2$	$36 \pm2$
	08$^\mathrm{h}$32$^\mathrm{m}$08$^\mathrm{s}$.65	$-43^\circ13'49.''5$	$6.9 \pm0.4 \times3.3 \pm0.2$	$111 \pm3$	$0.95 \pm0.05$	$196 \pm11$

Molecular Outflows

Figure 3(a) displays the outflow image of CO ($J = 2$–1) line emission. At least three distinct bipolar outflows are observed, originating from the NE and SW2 condensations, as well as the central protobinary system (C1 and C2). The SW1 and E condensations do not exhibit clear outflow features. However, weak blueshifted CO emission near SW1 may indicate a tentative outflow. The redshifted counterpart of the tentative SW1 outflow is difficult to identify due to blending with nearby redshifted outflows, and any detection remains highly uncertain. For the central bipolar outflow, it is challenging to identify the driving source owing to the close separation between C1 and C2. However, it is more likely to originate from C2, the primary component, as the continuum position angle of C2, $147 \pm1^\circ$, is nearly perpendicular to the outflow direction. In this case, the weaker redshifted outflow lobe to the northwest of C1 may represent the outflow from the C1 condensation; its blueshifted counterpart would then overlap with the blueshifted lobe of C2’s outflow. If so, this NW–SE outflow is nearly perpendicular to the continuum PA of C1, $60 \pm7^\circ$.

Figure 3(b) shows the outflow image of SiO ($J = 5$–4) line emission. Similar to CO, SiO is also a tracer of molecular outflows from MYSOs ((M. Liu et al. 2021; N. Guerra-Varas et al. 2023)). SiO often traces central jet-like features along the axes of broader CO outflows (e.g., (J. A. López-Vázquez et al. 2024)). Consistent with this expectation, the detected blueshifted SiO emission ridge is located near the center of the CO lobe, and its extension intersects with condensation C2. This supports the earlier interpretation that the central outflow originates from the primary component, C2. There is also redshifted SiO emission corresponding to the redshifted CO lobe northwest of C1. However, its blueshifted counterpart is blended with the outflow emission from C2; therefore, higher angular resolution observations are required to disentangle these components and confirm the presence of an outflow originating from C1.

We derived the physical and dynamical parameters of the three outflows using a similar approach to that of H. Beuther et al. (2002) (see also (S. Cabrit & C. Bertout 1990; K.-T. Kim & S. E. Kurtz 2006)). We assumed that the excitation temperature of CO gas is 30 K (H. Beuther et al. (2002)), and the line intensity ratio of ${}^{12}\mathrm{CO}$ to ${}^{13}\mathrm{CO}$ in the line wing velocity range is 10 (M. Choi et al. (1993)). The ${}^{13}\mathrm{CO}$ ($J = 2$–1) line was also covered in the QUARKS frequency band, but for this target, the line emission is weak in the high-velocity wings, making it difficult to calculate a reliable line intensity ratio. Therefore, we adopted a typical literature value instead. The assumed fractional abundance of $^{13}\mathrm{CO}$ is $[\mathrm{H}_2]/⁷ = 8.9 \times 10^{5} (S. Cabrit & C. Bertout 1992), and the mean molecular weight is 2.8. For simplicity, we took the statistical mean value of the inclination angle (i), assuming that the outflows are randomly oriented. For a random distribution of orientations, the statistical mean inclination angle is not 45^\circ$but$\sim57\mathord{.}3^\circ (S. Bontemps et al. 1996), because the probability distribution is weighted by \sin i, giving greater weight to nearly edge-on orientations. We first derived the maximum three-dimensional velocity of each outflow lobe, v_{\max}$, defined as$(\lvert v_{\mathrm{LSR}} - v_{\mathrm{sys}}\rvert/\cos i)_{\max}$, where$v_{\mathrm{LSR}} is the observed LSR velocity (see Table 4). Dynamical parameters, including momentum, energy, and mechanical force, were calculated following the v_{\max} method described in Section 4.2.1 of N. van der Marel et al. (2013). This approach assumes that each outflow lobe moves at a single characteristic bulk velocity v_{\max}. Although this method may underestimate the outflow force by a factor of 1.5–3 compared to alternative techniques discussed in N. van der Marel et al. (2013), it remains one of the most widely adopted approximations in the literature (S. Cabrit & C. Bertout 1992; M. R. Hogerheijde et al. 1998; T. P. Downes & S. Cabrit 2007). H. M. de Villiers et al. (2014) summarize in their Table 4 how these dynamical parameters are affected by the inclination angle. The accretion rate was estimated based on the conservation of momentum between the jet and the outflow. Due to the non-detection of a high-velocity SiO component (Figure 3), which precluded a direct measurement of the jet velocity, we adopted a representative value of v_{\mathrm{jet}} = 500\ \mathrm{km\ s^{-1}} (J. Martí et al. 1998). This approach differs from that of H. Beuther et al. (2002), which assumed a fixed velocity ratio of v_{\mathrm{jet}}/v_{\mathrm{out}} = 20$based on the averaged outflow velocity ($v_{\mathrm{out}} \sim 28\ \mathrm{km\ s^{-1}}) derived from their single-dish data. Since our ALMA observations resolve much higher outflow velocities (v_{\mathrm{out}} > 100\ \mathrm{km\ s^{-1}}), adopting this constant ratio would result in an unphysically large v_{\mathrm{jet}}$, leading to an underestimation of$\dot{M}_{\mathrm{jet}}$and$\dot{M}_{\mathrm{acc}}$. We therefore employed a fixed$v_{\mathrm{jet}} consistent with typical high-mass protostellar jets (J. Martí et al. 1998). Additionally, we assumed a ratio of jet mass-loss rate to accretion rate (\dot{M}_{\mathrm{jet}}/\dot{M}_{\mathrm{acc}}) of approximately 0.3 (K. Tomisaka 1998; F. H. Shu et al. 1999). Due to the various assumptions, the resulting outflow parameters have large uncertainties: a factor of \sim3$for the outflow mass, a factor of$\sim10$for the mechanical force, and a factor of$\sim30 for the mechanical luminosity (H. M. de Villiers et al. 2014). Table 4 lists the derived parameters. Our estimates are consistent with those of MYSOs with \dot{M}_{\mathrm{out}} \ge 10^{-4}\ M_\odot\ \mathrm{yr}^{-1} (A. Y. Yang et al. 2018; T. Baug et al. 2021), but they are significantly larger than those of low-mass YSOs with \dot{M}_{\mathrm{out}} \le 10^{-6}\ M_\odot\ \mathrm{yr}^{-1}$ (M. M. Vazzano et al. 2021; S. Dutta et al. 2024), indicating that the driving sources are most likely MYSOs.

	$v_{\mathrm{LSR,b}}^{a}$	$v_{\mathrm{LSR,r}}^{a}$	$N_b$	$N_r$	$M_b$	$M_r$	$M_{\mathrm{out}}$	$p$	$E$	$r$	$t$	$\dot{M}_{\mathrm{out}}$	$\dot{M}_{\mathrm{acc}}$	$F_m$	$L_m$
NE	[-56.2, 9.2]	[18.7, 62.5]	3.02	3.23	4.88	2.03	6.91	8.15	9.77	8.50	367	1.88	1.38	2.22	2.21
C	[-68.3, 9.2]	[18.7, 55.6]	5.25	6.16	19.20	10.63	29.83	37.37	50.31	9.39	388	7.68	5.87	9.62	10.7
SW2	[-59.4, 9.2]	[18.7, 49.2]	3.34	4.67	7.42	31.46	38.87	29.68	26.17	19.90	939	4.14	2.77	3.16	2.31

Circumbinary Disk Candidate around the Protobinary System

The $\mathrm{CH}_3\mathrm{CN}$ molecular lines have frequently been used to observe accretion disks around MYSOs (e.g., K. G. Johnston et al. 2015; J. D. Ilee et al. 2018). The line emission was detected toward all the six continuum condensations, although it was only marginally detected in the E condensation (Table 2). Figure 4 shows the integrated line intensity maps of several molecular lines, including $\mathrm{CH}_3\mathrm{CN}$, HNCO, and $\mathrm{SO}_2$, which trace an elongated structure surrounding the C1 and C2 condensations. J. D. Ilee et al. (2016) found that molecules such as $\mathrm{CH}_3\mathrm{OH}$ and $\mathrm{H}_2\mathrm{CO}$ also trace the circumstellar disk in the case of the MYSO G11.92−0.61 MM1, similar to $\mathrm{CH}_3\mathrm{CN}$. Figure 5 suggests that this elongated structure may represent a circumbinary disk, as a consistent velocity gradient is observed along its major axis in multiple molecular lines. This finding is similar to the ALMA observational results of IRAS 16547−4247 in which $\mathrm{CH}_3\mathrm{CN}$ traces the circumbinary disk of a high-mass protobinary system (K. E. I. Tanaka et al. 2020). Also, the higher-velocity component relative to the systemic velocity is located near the center, consistent with Keplerian rotation. Additionally, enhanced $\mathrm{CH}_3\mathrm{CN}$ velocity dispersion is largely confined to the extent of the proposed circumbinary disk, apart from a localized enhancement immediately west of the C1/C2 center that may indicate a minor outflow contribution (see Appendix § 8). Assuming that the circumbinary disk is circular, we derived its inclination angle from the six transitions shown in Figure 4 (a)–(Figure 4f) as $\cos^{-1}(D_{\mathrm{min}}/D_{\mathrm{maj}}) = 55 \pm3^\circ$, where $D_{\mathrm{min}}$ and $D_{\mathrm{maj}}$ are the minor and major diameters, respectively. The small uncertainty displays that the six measurements agree well. Interestingly, the majority of the detected molecular lines, with the exception of the $\mathrm{SO}_2$ lines, are significantly brighter in C1 than in C2, which is opposite to what is observed in the continuum emission. Another interesting aspect of our findings is that $\mathrm{CH}_3\mathrm{OH}$, $\mathrm{H}_2\mathrm{CO}$, $\mathrm{CH}_3\mathrm{CN}$, and HNCO exhibit an arc-like structure extending southeast from C1. This feature may represent a streamer, similar to those observed in the high-mass protostellar core G336.01−0.82 (F. A. Olguin et al. 2023, 2025). We discuss this feature further in Section 4.3.

Disk Candidate Associated with the Other Condensations

The central protobinary system, C1 and C2, is not the only outflow-driving source in IRAS 08303–4303; the NE and SW2 condensations also drive molecular outflows. Since bipolar outflows are closely associated with accretion disks, we might expect to detect disks in these components as well. The NE condensation does not show rich molecular line features, while the SW2 condensation displays several COM lines, although the lines are much weaker than those observed in the central region (see Figure 2 and Table 2). Figure 6 demonstrates the velocity structure of the CH$_3$OH line around the SW2 condensation. There is a velocity gradient in the northeast–southwest direction, which is roughly parallel to the major axis of SW2 but perpendicular to the associated CO bipolar outflow. Therefore, considering that the SW2 condensation is about 1600 × 500 au in size (Table 3), it seems to be a circumstellar disk around the embedded high-mass (proto)star. Since this condensation is marginally resolved in the QUARKS observations, higher angular resolution (and sensitivity) observations are required to investigate its velocity structure in detail.

Discussion

Kinematics of the Circumbinary Disk

The total binary mass (enclosed mass, $M_{\mathrm{enc}}$) can be derived from the rotation kinematics of the circumbinary disk combined with the inclination information. The inclination of the circumbinary disk was determined earlier by fitting the integrated line intensity map, as shown in Figure 4, yielding 55°. Assuming Keplerian rotation, the enclosed mass is given by

where $v$ is rotation velocity at radius $r$, $G$ the gravitational constant. The systemic velocity of the central protobinary system was measured to be 13.7 km s$^{-1}$ by fitting the 234.187 GHz SO$_2$ line profile, which shows similar intensities toward the C1 and C2 condensations (see Figure 4). For other molecules exhibiting asymmetric intensity distributions, the fitted peak velocities are considerably closer to the velocity of C1 than to the actual systemic velocity. The rotation center is defined as the point with the systemic velocity in the intensity-weighted velocity map along the line connecting the C1 and C2 continuum peaks. Using this information, together with the CH$_3$CN ($J = 12$–11 $K = 3$) position–velocity (PV) diagram in Figure 7, we employed a Keplerian velocity distribution model fitter (Keplರ್fit) to derive $M_{\mathrm{enc}}$ ((F. Bosco et al. 2019, 2023)). A 7σ emission edge was used for the fitting after excluding the inner region within [−1000, +1000] au. Because the circumbinary disk does not contain actual mass at the rotation center, the rotation at small radius may deviate from Keplerian motion, and this region was therefore masked. The outer region beyond 3000 au was also masked because of the limited velocity resolution (∼ 1.3 km s$^{-1}$), which prevents the velocity structure at low velocities from being resolved. The data points used in Keplerfit are marked with black crosses in Figure 7(b). The inclination, systemic velocity, and rotation center were fixed, while the mass was set as the only free parameter. The enclosed mass was estimated to be 61 ± 5 $M_{\odot}$ from the best fit, which is delineated by the green dashed lines.

We also performed an additional fitting of the CH$_3$CN ($J = 12$–11 $K = 3$) PV diagram using the method of (G. M. Williams et al. 2022). This fitting considers both Keplerian rotation along the tangential direction and free-fall motion in the radial direction. The LSR velocity within the disk is then expressed as

where $x$ and $z$ are coordinates centered at the rotation center. Here, $x$ lies in the disk plane, perpendicular to the line of sight, while $z$ also lies in the disk plane but perpendicular to $x$. The distance from the rotation center is $R=\sqrt{x^2+z^2}$, which is constrained to the range $[R_{\mathrm{in}},R_{\mathrm{out}}]$. We fixed $v_{\mathrm{sys}}=13.7\ \mathrm{km\ s^{-1}}$, $i=55^\circ$, and $[R_{\mathrm{in}},R_{\mathrm{out}}]=[1000,3000]\ \mathrm{au}$, consistent with the Keplerfit procedure, and treated the mass as the only free parameter. We then fit $v_{\mathrm{LSR}}$ as a function of $x$ using the white dots in Figure 7(b), selected from the $7\sigma$ emission edge within 3000 au. We found that a model with an enclosed mass of $32\pm2\ M_\odot$ best reproduces the data. The fitting result is shown as the red solid line in Figure 7(b). Having determined that the enclosed mass range is between 32 and 61 $M_\odot$ according to the two methods, we can derive the individual component masses by measuring the distance from the rotation center to each condensation peak. The C1 peak is located 630 au from the rotation center, while the C2 is 270 au away. By definition of the center of mass, the distance ratio from the rotation center is the inverse of the mass ratio. Based on the observed distance ratio of 2.3:1, we estimated a mass ratio of 1:2.3. This yields a primary star mass in C2 of 22–43 $M_\odot$ and a secondary star mass in C1 of 10–18 $M_\odot$.

We assessed the gravitational stability of the circumbinary disk using the central stellar mass, 1.3 mm continuum emission, and temperature information, specifically through the Toomre $Q$ parameter (A. Toomre 1964). This parameter evaluates the balance between the local gravitational forces of the disk and the stabilizing effects of shear and thermal pressure [2]. Following A. Ahmadi et al. [2], we write $Q$ as

Here, $c_s$ denotes the sound speed within the disk, $\Omega$ is the Keplerian angular velocity, and $\Sigma$ is the surface density of the disk. The adiabatic index is $\gamma=7/5$, $k_{\mathrm{B}}$ is the Boltzmann constant, $\mu$ is the mean molecular weight of the gas (assumed to be 2.8), and $m_{\mathrm{H}}$ is the mass of a hydrogen atom. $M_\ast$ represents the mass of the central protostar(s), while $M_{\mathrm{disk}}(\leq r)$ is the enclosed disk mass within radius $r$, which is neglected here because the stellar-to-disk mass ratio is very high ($M_\ast/M_{\mathrm{disk}}\gtrsim15$). Finally, $S_\nu$ is the continuum intensity, and $\theta_{\mathrm{B}}$ is the beam solid angle. If the Toomre $Q$ parameter is smaller than 2, the disk is gravitationally unstable. It is also roughly proportional to $M_\ast/M_{\mathrm{disk}}$, and disks with masses larger than 10–20% of the stellar mass tend to be gravitationally unstable (K. Kratter & G. Lodato 2016; A. Ahmadi et al. 2023). So, with $M_{\mathrm{enc}}\sim32$–61 $M_\odot$ and a relatively small circumbinary disk mass $M_{\mathrm{disk}}\sim2\ M_\odot$, this circumbinary disk is overall gravitationally stable. However, Figure 7(c), the distribution of the calculated Toomre $Q$ parameter for $M_{\mathrm{enc}}=32\ M_\odot$, indicates that certain regions—especially southeast of C2—are gravitationally unstable ($Q<2$). Additionally, given that the disk midplane may be significantly colder than the gas temperature used in the calculation, other regions could also be gravitationally unstable. For $M_{\mathrm{enc}}=61\ M_\odot$, the values would be approximately 1.5 times higher.

Comparison with Previous Studies

To date, only a few high-mass protobinary systems have been identified. IRAS 20126+4104 was the first of these systems, consisting of a primary (proto)star with a mass of $\sim7\ M_\odot$ and a companion located at a projected separation of around 850 au (T. K. Sridharan et al. 2005). Another candidate high-mass binary system is G11.92–0.61 MM1 with a projected separation of about 2000 au (J. D. Ilee et al. 2018). The primary component has a mass of roughly $\sim34\ M_\odot$, while the secondary is a low-mass companion ($<0.6\ M_\odot$), likely formed through disk fragmentation (G35.20–0.74 N–B contains a binary system with a companion mass of $\sim6\ M_\odot$, an apparent separation of roughly 800 au from the primary component with a mass $\gtrsim11\ M_\odot$, and is surrounded by a circumbinary disk traced by CH$_3$CN (Á. Sánchez-Monge et al. 2013; M. T. Beltrán et al. 2016)). The IRAS 16547–4247–E protobinary system has a total enclosed mass of about 25 $M_\odot$ and an apparent separation of $\sim300$ au (L. A. Zapata et al. 2019). For this system, individual circumstellar disks as well as the circumbinary disk have been identified (K. E. I. Tanaka et al. 2020). W3(H$_2$O) is another example of a high-mass protobinary system where both a circumbinary disk and individual circumstellar disks have been observed [1]. This system consists of two $\sim10\ M_\odot$ protostars at a projected separation of $\sim2300$ au, surrounded by a circumbinary disk spanning $\sim6000$ au in diameter. Notably, the kinematics of its circumbinary disk deviates significantly from typical Keplerian rotation. NGC 7538 IRS1 is a high-mass binary system with a projected separation of $\sim430$ au, and its total luminosity corresponds to a stellar mass of about 30 $M_\odot$. (H. Beuther et al. 2017) IRAS 07299–1651 features a close binary system with component masses of 12.5 $M_\odot$ and 10 $M_\odot$, with an apparent separation of only ∼180 au (Y. Zhang et al. 2019). Finally, IRAS 17216–3801 is a well-characterized high-mass binary system with component masses of approximately 20 $M_\odot$ and 18 $M_\odot$, with a circumbinary disk and a close projected separation of ∼170 au, as observed in near-infrared emission (S. Kraus et al. 2017).

Among the aforementioned high-mass protobinary systems, IRAS 16547–4247–E is the most extensively studied target regarding its circumbinary disk. L. A. Zapata et al. (2019) identified a circumbinary disk traced by CH$_3$OCHO and determined its enclosed mass through kinematic analysis to be 25 ± 3 $M_\odot$ under the assumption of Keplerian rotation at an inclination of 55 ± 5°. The radius of this disk is 1190 ± 290 au at an adopted distance of 2.9 kpc. Higher-resolution (∼0.′′05) observations by K. E. I. Tanaka et al. (2020) revealed two individual circumstellar disks in the H$_2$O and refractory molecular (e.g., SiS, NaCl, and SiO isotopologues) lines, which are rotating roughly in the opposite directions. Besides, the circumbinary disk was observed by traditional disk tracers such as CH$_3$CN and SO$_2$ similarly to IRAS 08303–4303–C in this study (Figure 5). Unfortunately, we cannot detect any circumstellar disk within the circumbinary disk of IRAS 08303–4303–C mainly due to the limited angular resolution of the QUARKS observations.

The protobinary system of IRAS 08303–4303–C has a larger separation of ∼900 au, a larger disk radius of about 3000 au, and a larger enclosed mass of 32–61 $M_\odot$, compared to the IRAS 16547–4247–E. Thus, it can be considered a scaled-up version of IRAS 16547–4247–E. However, despite its larger enclosed mass, the bolometric luminosity of IRAS 08303–4303 (∼10$^4$ $L_\odot$) is one order of magnitude smaller, compared to the ∼10$^5$ $L_\odot$ of IRAS 16547–4247.

If a primary component (C2) has a mass between 22 and 43 $M_\odot$, its bolometric luminosity is expected to be similar to or greater than 10$^5$ $L_\odot$ at the zero-age main sequence (ZAMS) (W. D. Vacca et al. 1996; F. Martins et al. 2005; B. Davies et al. 2011). However, its actual observed luminosity (∼7000 $L_\odot$) is at least an order of magnitude lower. One explanation for this discrepancy is that a high accretion rate has caused an increase in the stellar radius while decreasing the surface temperature and luminosity. According to T. Hosokawa & K. Omukai (2009) and R. Kuiper & H. W. Yorke (2013), the “bloated” phase of massive protostars is characterized by a high accretion rate, ranging from a few 10$^{-4}$ to a few 10$^{-3}$ $M_\odot$ yr$^{-1}$, alongside a stellar radius of several 100 $R_\odot$. From § 3.2, the mass accretion rate calculated from CO outflows exceeds 5 × 10$^{-4}$ $M_\odot$ yr$^{-1}$, and could be even higher if the outflow inclination angle is greater than 57°, which is consistent with this “bloating” phase. G. M. Williams et al. (2022) also discussed this phenomenon in relation to the high stellar mass ($40$–$70\,M_\odot$) observed in G19.01–0.03 MM1, which has a moderate luminosity of approximately $10^4\,L_\odot$, a situation similar to that presented in this work. In addition, the presence of unresolved components could contribute to this discrepancy. For instance, according to Table 1 in B. Davies et al. (2011), a single $30\,M_\odot$ star has $L \gtrsim10^5\,L_\odot$, whereas two $15\,M_\odot$ stars have $L < 4 \times10^4\,L_\odot$. If the estimated primary protostellar mass in C2 is distributed among two or more unresolved sources, the resulting luminosity would be lower, consistent with the observations. In fact, this system is further resolved into multiple individual components in recent higher-resolution observations K. Yang et al. (2026), which favors this possibility although it is not known how many of the components contain embedded (proto)stars.

Molecular Streamer

Recently, nonspherical inflows referred to as streamers have been discovered both in low- and high-mass star-forming regions. These streamers transport materials from the inner envelope to the protostellar disk (M. T. Valdivia-Mena et al. 2022; M. Fernández-López et al. 2023; X.-C. Liu et al. 2025). We also detected a streamer candidate around the circumbinary disk. Figure 4 illustrates an arc-like structure extending from the C1 condensation, which is observed in various molecules such as CH$_3$CN, HNCO, CH$_3$OH, and H$_2$CO. The H$_2$CO has been utilized to trace streamers around low-mass YSOs (M. T. Valdivia-Mena et al. 2022), while CH$_3$OH and CH$_3$CN have been employed to trace streamers around MYSOs (M. Fernández-López et al. 2023; F. A. Olguin et al. 2023, 2025). A velocity gradient of $0.6$ km s$^{-1}$ ($10^3$ au)$^{-1}$/$\sin i$ is observed along the streamer, accompanied by an increase in velocity dispersion toward C1 (Figure 8). A binary system associated with a streamer resembles the low-mass binary system [BHB2007] 11, where spiral-like infalling streams link the circumbinary disk and the envelope (F. O. Alves et al. 2017). However, such smaller-scale streamers, which connect the circumbinary disk to the individual circumstellar disks in [BHB2007] 11 (F. O. Alves et al. 2019), are not detected here presumably due to the coarse resolution. The low-mass triple system IRAS 04239+2436 also exhibits spiral streamers produced by gravitational interactions among its members (J.-E. Lee et al. 2023).

The C1 condensation is weaker than C2 in 1.3 mm continuum emission, while it exhibits higher integrated intensities than the latter in most molecular lines, as in Figure 2. This difference may be attributed to differences in temperature and/or molecular abundance. To investigate this further, we made a rotational temperature map using the intensities of multiple CH$_3$OH transitions, employing the assumption of optically thin line emission in LTE (for detailed methodology, see Appendix § 7). The results indicate that the region surrounding C1 has a higher temperature compared to C2 (Figure 8). The temperature of the streamer ranges from 50 to 150 K, increasing as it approaches C1. This value is much higher than the 10–15 K typically found in low-mass streamers (J. E. Pineda et al. 2020; M. Kido et al. 2025), and comparable to the ∼100 K high-mass streamer reported by F. A. Olguin et al. (2025). The molecular abundances toward C1 are notably higher than those of C2; specifically, the abundances of CH$_3$CN and CH$_3$OH are 5–7 times greater (see Appendix § 7). One possible explanation for C1’s elevated temperature and molecular richness is the impact of materials accreting via the streamer onto the disk plane, which may lead to accretion shock.

It would be worthwhile to calculate both the mass of the streamer and the mass accretion rate through it. However, the spiral structure resembling the streamer is not clearly visible in the 1.3 mm continuum emission. This suggests that the weak continuum emission detected in the streamer region may originate from an envelope along the line of sight, rather than being directly associated with the streamer itself. Consequently, calculating the mass from continuum observations is challenging. Similarly, estimating the mass using CH$_3$OH or H$_2$CO—molecular lines that trace the streamer—is also difficult, as their abundances vary strongly with environmental conditions. Another issue is that performing an accurate kinematic analysis is complicated by the unknown inclination of the streamer. If the streamer lies in the plane of the disk, we could assume an inclination of $55^\circ$, but this assumption is problematic, as there also been reports of streamers located off the disk midplane (X.-C. Liu et al. 2025). Future higher-sensitivity observations would facilitate mass estimation, while higher angular- and velocity-resolution studies would enable a detailed analysis of the streamer and its interface with the disk.

Nature of the IRAS 08303–4303

Binary systems generally form through two primary mechanisms: disk fragmentation and turbulent fragmentation. In the former, gravitational instability within a disk leads to the formation of close binaries ($<500$ au), typically resulting in aligned circumstellar disks (I. A. Bonnell & M. R. Bate 1994; K. M. Kratter et al. 2010; F. O. Alves et al. 2019). In contrast, wide binaries (> 500 au) usually form through turbulent fragmentation, which produces systems with misaligned circumstellar disks due to the initially random angular momentum distribution in a turbulent core (S. S. R. Offner et al. 2010; J.-E. Lee et al. 2017). For instance, M. R. Bate (2018) showed in clump-scale numerical simulations that multiple systems formed in such environments often possess misaligned disks, some of which are associated with circumbinary material. These disordered systems are consistent with observations of the binary protostar IRS 43, which comprises a circumbinary disk encompassing two misaligned circumstellar disks ($\gtrsim60^\circ$), attributed to either turbulent fragmentation or recent dynamical perturbations (C. Brinch et al. 2016). Despite such initial misalignment, a circumbinary disk can form a posteriori through the subsequent accretion of gas from the surrounding envelope. M. R. Bate (2018) also demonstrated that the orientation of a circumbinary disk is governed by the angular momentum of the infalling material, even when it is significantly misaligned with the existing binary systems. Circumbinary disks form naturally around binary systems when the specific angular momentum of the infalling material is sufficiently large (M. R. Bate & I. A. Bonnell 1997; T. Matsumoto et al. 2019). In such cases, circumbinary disks are established on short timescales, typically within a few orbital periods of the binary system (Y. Ochi et al. 2005; T. Matsumoto et al. 2019).

The projected separation of $\sim900$ au in IRAS 08303–4303–C indicates that this system is a wide binary. The PAs of the C1 and C2 condensations are nearly perpendicular (Table 3). While the PA of an individual continuum condensation may not perfectly reflect the true orientation of the circumstellar disk due to potential influence from unresolved substructures, the fact that the outflows from C1 and C2 also appear to be oriented in approximately perpendicular directions further corroborates this misalignment. Such misaligned disk–outflow systems suggest that this wide binary formed through turbulent fragmentation on a large scale. The circumbinary disk might form subsequently after the binary formation. On the other hand, additional unresolved components could have formed within the circumbinary disk due to disk fragmentation. Notably, the C2 condensation appears highly elongated toward the southeast from its continuum peak (see, e.g., Figure 7(a)). This non-Gaussian feature hints the presence of an unresolved component on the southeastern side. Moreover, the low Toomre $Q$ parameter in this region (see Figure 7(c)), along with the discrepancies in the observed luminosity and the derived enclosed mass discussed in § 4.2, also suggests the possible existence of an unresolved component.

As noted above, K. Yang et al. (2026) resolved the C1 and C2 condensations into at least three distinct components each, one of which lies southeast of the C2 peak. While the limited resolution of our QUARKS data precluded detailed analysis of individual disk kinematics, they marginally detected a circumstellar disk candidate toward one component within C1 in $\mathrm{H_2O}$ line emission, which often serves as a more reliable tracer for hot, inner disk dynamics than COMs (T. Hirota et al. 2017; A. Ginsburg et al. 2018; K. E. I. Tanaka et al. 2020. However, their molecular line analysis focused solely on the single component. Therefore, a joint analysis of our comprehensive data with higher-resolution data, including theirs, is needed to characterize disk rotation in the other components and accretion flows onto individual circumstellar disks. Although the continuum condensation and outflow morphologies observed in this study appear to support the misalignment among the individual components, unresolved substructures might produce some uncertainty in the measured orientations of the innermost disks and their associated outflows. Such a joint analysis would help measure their true orientations and test our interpretation of the misalignment and turbulent fragmentation.

Summary

We carried out ALMA observations of the high-mass star-forming region IRAS 08303–4303 in Band 6 with an angular resolution of about 0\farcs3, as part of the QUARKS project. We detected multiple dense condensations and the associated disk-outflow systems, and investigated the physical and dynamical properties of them. The main results are summarized as follows.

1. We detected six dense continuum condensations at 1.3 mm. At least four of these condensations (NE, C1, C2, SW2) show signs of star formation activity, as indicated by the line-rich spectra and associated outflows. The NE, C2, and SW2 condensations demonstrate clear bipolar molecular outflows, while the C1 (and possibly SW1) condensation is associated with weak outflow features.

2. The C1 and C2 condensations, located apparently $\sim900$ au apart, form a wide protobinary system and reside within a single rotating structure. This structure is characterized as a ‘circumbinary disk’ and exhibits a velocity pattern similar to Keplerian rotation, with higher velocities observed near the center. Multiple molecules, including CH$_3$CN, SO$_2$, CH$_3$OH, HC$_3$N, HNCO, and OCS, trace this rotation.

3. The total mass of the protobinary system was derived using both Keplerian rotation and combined Keplerian rotation–infall models. The total mass estimates range between 32 and 61 $M_\odot$, according to two different fitting results. The masses of the primary and secondary protostars were computed to be 22–43 and 10–18 $M_\odot$, respectively, based on their distances from the rotation center.

4. The temperature of the circumbinary disk was determined from a rotational diagram analysis using multiple transitions of CH$_3$OH, yielding a value of approximately 130 K. The mass of the circumbinary disk was then estimated to be about 2 $M_\odot$ based on the continuum flux and temperature data. Given the relatively low disk mass compared to the high protostellar mass, the disk is overall gravitationally stable, except for the southeastern region of the C2 condensation showing $Q < 2$.

5. From the outflow parameters, the mass accretion rate onto the primary protostar in C2 was derived to be greater than $5 \times10^{-4}$ $M_\odot$ yr$^{-1}$. This high accretion rate may explain the discrepancy between the observed bolometric luminosity of less than $10^4$ $L_\odot$ and the relatively high estimated protostellar mass of 22–43 $M_\odot$. During the protostellar stage, a high accretion rate can cause the stellar surface to expand, resulting in a lower luminosity compared to the zero-age main sequence. In addition, the possible presence of an unresolved component around C2 could also account for this inconsistency.

6. The C1 condensation is brighter in most molecular lines and exhibits a higher temperature compared to the primary component, C2. This may be attributed to a streamer connected to C1, which is observed in several molecular lines, including CH$_3$OH, H$_2$CO, CH$_3$CN, and HNCO. Accreting material flowing along the streamer may interact with the disk plane, leading to enhanced temperatures and molecular abundances around C1 due to accretion shock.

The IRAS 08303–4303 region provides valuable insights into the rarely observed circumbinary disks around massive YSOs. It serves as an excellent laboratory for future studies aimed at understanding high-mass multiple formation through higher-resolution observations.

Acknowledgments

We thank the anonymous referee for his/her constructive and thoughtful comments. This paper makes use of the following ALMA data: ADS/JAO. ALMA#2019.1.00685.S and 2021.1.00095.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan, China), and KASI(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. GG and LB gratefully acknowledge support by the ANID BASAL project FB210003. HLL is supported by the National Natural Science Foundation of China (No. 12673033), by Yunnan Fundamental Research Project (grant No. 202301AT070118, 202401AS070121), and by Xingdian Talent Support Plan-Youth Project. DL is supported by the National Natural Science Foundation of China (NSFC) Programs No. 12588202. SRD acknowledges support from the Fondecyt Postdoctoral fellowship (project code 3220162) and ANID BASAL project FB210003. Suinan Zhang acknowledges support from the Strategic Priority Research Program of the Chinese Academy of Sciences (CAS) Grant No. XDB0800303 and the NAOJ ALMA Scientific Research Grant Code 2025-29B. MJ acknowledges the support of the Research Council of Finland Grant No. 348342.

Facilities: ALMA

Software: CASA (CASA Team et al. 2022), Keplerfit (F. Bosco et al. 2019, 2023)

Temperature calculation

Assuming the LTE condition, the rotational temperature ($T_{\mathrm{rot}}$) can be calculated using the equations provided by P. F. Goldsmith & W. D. Langer (1999), as shown below:

Here $\nu$ is the rest frequency; $W$ the integrated line intensity; $h$ the Planck constant; $c$ the speed of light; $A_{ul}$ the Einstein coefficient; $N_{u}^{\mathrm{thin}}$ the upper-level column density when the lines are optically thin; $g_{u}$ the upper-level degeneracy; $C_{\tau}=\tau/(1-e^{-\tau})$ the optical depth correction factor ($C_{\tau}\simeq1$ when $\tau\ll1$); $N$ the total column density of the molecule; $Z$ the partition function; and $E_{u}$ the upper-level energy. From Table 2, nine transitions of CH$_3$CN and six of CH$_3$OH were detected within the QUARKS frequency band, each corresponding to a different upper-level energy. By plotting $\ln(N_{u}^{\mathrm{thin}}/g_{u})$ against $E_{u}/k_{\mathrm{B}}$, the inverse of the slope yields the rotation temperature, while the y-intercept allows us to determine the total column densities of CH$_3$CN and CH$_3$OH (see Figure A1). For CH$_3$OH, we performed the rotation diagram analysis by combining the transitions of A-type and E-type species. Similar methodologies have been adopted in several recent studies where the number of detected CH$_3$OH lines was insufficient to determine the A/E ratio independently (E. Bianchi et al. 2020; G. Baek et al. 2022; P. Freeman et al. 2023). This simplification is further justified by the high-temperature environment ($T>37$ K) where the A/E ratio is expected to approach unity within 0.5% (E. S. Wirström et al. 2011). The hot-core nature of our target ($T>100$ K) significantly exceeds the A/E equilibration threshold, supporting the single excitation temperature assumption for both species. Molecular spectroscopic data (e.g., rest frequencies, Einstein coefficients, upper-level degeneracies, and partition functions) were obtained from the CDMS catalog (H. S. P. Müller et al. 2001, 2005; C. P. Endres et al. 2016). Partition functions were interpolated to the appropriate temperatures.

All the rotational temperature and molecular column density calculations for C1, C2, SW2, and the overall central protobinary system were performed based on the optically thin assumption, where $C_{\tau}$ is considered to be unity for all lines. The CH$_3$CN $J=12$–$11$ $K=0$ and $K=1$ transitions are partly blended due to their close rest frequencies, and the CH$_3$CN $J=12$–$11$ $K=6$ transition overlaps with nearby CH$_3^{13}$CN and HNCO lines. Therefore, Gaussian fitting was used to separate the integrated line intensities of these transitions. For the SW2 condensation, the calculations were only conducted using CH$_3$OH lines because of the weak intensities of the CH$_3$CN lines. Assuming that the dust temperature ($T_{\mathrm{d}}$) is equal to the derived rotational temperature, the average H$_2$ column density within the source solid angle $\theta$ was estimated from the 1.3 mm continuum flux density $F_{\nu}$ using the following equation from J. Kauffmann et al. (2008):

Molecular abundance was calculated as the ratio of the molecular column density to the H$_2$ column density. From the temperatures and molecular column densities derived from the rotation diagram analysis, we also estimated the optical depths of several CH$_3$OH and CH$_3$CN transitions using the non-LTE radiative transfer code RADEX (F. F. S. van der Tak et al. 2007). The CH$_3$OH line emission is more optically thin than that of CH$_3$CN toward all regions; therefore, we present the temperature and mass estimates derived from CH$_3$OH in Table 3. All results from the rotation diagram analysis are summarized in Table A1. Figure A2 shows the maps of temperature, molecular column density, and molecular abundance obtained by applying equation (5) to each pixel of the integrated line intensity maps of CH$_3$OH, where at least four CH$_3$OH transitions were detected above $5\sigma$.

Condensation	$T_{\rm rot}$ (K)	$N_{\rm CH_3CN}$ ($10^{15}$ cm$^{-2}$)	$N_{\rm H_2}$ ($10^{23}$ cm$^{-2}$)	$X_{\rm CH_3CN}$ ($10^{-8}$)	$\tau_{J=12-11\ K=3}$	$\tau_{J=12-11\ K=5}$
C1	$297 \pm8$	$10.0 \pm0.6$	$2.7 \pm1.0$	$3.8 \pm1.4$	$1.43$	$0.44$
C2	$201 \pm7$	$3.4 \pm0.3$	$5.2 \pm1.9$	$0.6 \pm0.2$	$0.74$	$0.19$
SW2
Central	$175 \pm6^{a}$	$1.26 \pm0.09^{a}$	$1.1 \pm0.4$	$1.2 \pm0.4$	$0.33$	$0.08$
	CH$_3$OH
	$T_{\rm rot}$ (K)	$N_{\rm CH_3OH}$ ($10^{16}$ cm$^{-2}$)	$N_{\rm H_2}$ ($10^{23}$ cm$^{-2}$)	$X_{\rm CH_3OH}$ ($10^{-6}$)	$\tau_{42,3-31,2}$	$\tau_{102,9-93,6}$
C1	$235 \pm3$	$45 \pm1$	$3.4 \pm1.2$	$1.4 \pm0.5$	$0.79$	$0.38$
C2	$171 \pm6$	$14 \pm1$	$6.1 \pm2.2$	$0.22 \pm0.08$	$0.42$	$0.18$
SW2	$262 \pm9$	$9.1 \pm0.6$	$1.3 \pm0.5$	$0.7 \pm0.3$	$0.13$	$0.07$
Central	$127 \pm2$	$8.9 \pm0.3$	$1.5 \pm0.6$	$0.6 \pm0.2$	$0.48$	$0.16$

CH$_3$CN velocity dispersion

Figure B1 displays the velocity dispersion maps (Moment 2) of 3 CH$_3$CN transitions. Around the C2 condensation, all three transitions exhibit broad linewidths corresponding to FWHM $\gtrsim7$ km s$^{-1}$, assuming Gaussian line profiles. Unresolved circumstellar disk rotation may contribute to these broad linewidths. Toward the C1 condensation, the $K = 5$ transition shows a clearly larger velocity dispersion than the $K = 3$ and 4 transitions, which may indicate that it traces a hotter, faster-rotating region closer to the central source. Additionally, the $K = 3$ and 4 transitions exhibit more extended emission and reveal a localized enhancement in velocity dispersion immediately west of the C1/C2 center, which may be associated with the outflow. This possible outflow-related enhancement is not apparent in the more compact emission traced by the $K = 5$ transition.