Triple Spiral Arms of a Triple Protostar System Imaged in Molecular Lines

Although binary or multiple formation are the main modes of star formation (Offner et al. 2022), the main formation mechanism remains elusive. Observational and theoretical studies that zoom into young circumstellar and circumbinary disks in the early evolutionary stage are critical for understanding the formation process, before their history is erased through dynamical interactions (Lee et al. 2017; Moe & Kratter 2018). The interactions between the inner envelope and the central protostellar system, as well as the interactions among individual members, are important in multiple formation, since protostars grow their masses most aggressively in the early formation stage, when most of the mass still remains in the envelope. Therefore, simultaneous imaging of the inner envelope and the central compact multiple system is also required. One prominent structure that shows the interactions between forming binary and circumbinary material is the spiral arm structure (Takakuwa et al. 2017; Alves et al. 2019).

Two distinct formation mechanisms of multiple-star systems have been suggested: turbulent fragmentation for wide (>500 au) binaries (Offner et al. 2010; Lee et al. 2017) and disk fragmentation for close (<500 au) binaries (Kratter et al. 2010; Alves et al. 2019). However, at early stages, when the envelope is feeding material to the disk(s), the clean distinction breaks down. Instead, a hybrid process, where the turbulent envelope plays an important role in the disk fragmentation, the circumbinary/circumstellar disk structures, and the final masses of the individual central sources, may be important (Matsumoto et al. 2015b). In addition to turbulence, the magnetic field must also play an important role in the fragmentation process itself, and in the structures inside fragments, since cores with weak magnetic fields tend to form multiple systems (Machida et al. 2005).

To study the formation process of a close binary system when the environmental effect of the natal cloud may still remain in the system, we observed IRAS 04239+2436 with the Atacama Large Millimeter/submillimeter Array (ALMA) simultaneously in the 857 μm continuum and several molecular lines, with an angular resolution of ∼01 (∼14 au at 140 pc). IRAS 04239+2436 is an embedded close protobinary with a projected separation of 42 au at a distance of 140 pc (Reipurth et al. 2000). IRAS 04239+2436 has a luminosity similar to that of the Sun and very rich near-infrared (NIR) spectrum (Greene & Lada 1996). A high–spectral resolution NIR observation kinematically revealed a sub-astronomical-unit-scale gaseous disk with a Keplerian rotation around a central mass of 0.14 ∼ 0.2 M_⊙ (Lee et al. 2016). In the NIR images (2 μm) of IRAS 04239+2436 (Reipurth et al. 2000), the eastern source (Source A) is about 2.65 times brighter than the western source (Source B), and is thus considered as the primary component. In addition, a well-collimated northeast [Fe ii] 1.644 μm jet of HH 300 originates from Source A (Reipurth et al. 2000). Yet, the observed jet precesses at a short timescale (<30 yr), compared to the orbital period (∼270 yr) anticipated from the resolved binary, indicative of unresolved very close binary interactions (Reipurth et al. 2000). Therefore, IRAS 04239+2436 is very likely a triple system.

Here, we report ALMA observational results of a forming multiple protostellar system, IRAS 04239+2436 (Section 2), and its comparisons with numerical simulations (Section 3). We discuss the effect of the magnetic field on multiple formation and the implications of close multiple protostars for planet formation in Section 4, and the final conclusion is given in Section 5.

IRAS 04239+2436 was observed using ALMA during Cycle 3 (2015.1.00397.S; PI: Jeong-Eun Lee) on 2016 August 16 UT. Five spectral windows in Band 7 were set to cover several molecular lines, such as SO 8₈–7₇ (344.31061200 GHz), HCN 4–3 (354.50547590 GHz), HCO⁺ 4–3 (356.73422300 GHz), and SO₂ 10_4,6–10_3,7 (356.75518930 GHz). ¹⁰ The source velocity was V_LSR = 6.5 km s⁻¹ (Fuller & Ladd 2002). The bandwidths and spectral resolutions were 117.19 MHz and 122.070 kHz (∼0.1 km s⁻¹) for the SO line and 468.75 MHz and 244.141 kHz (∼0.2 km s⁻¹) for the other three. 40 12 m antennas were used in the C40-4 configuration, with baselines in the range from 21.4 m to 3.1 km. The total observing time on source was 29.8 minutes.

The data were initially calibrated using the CASA 4.7.0 pipeline (McMullin et al. 2007). The nearby quasar J0510+1800 was used as a bandpass and phase calibrator, and the quasar J0238+1636 was used as an amplitude calibrator. Self-calibration was applied for better imaging. The continuum image with Briggs weighting (robust = 0.5) is used for the self-calibration. Two phase calibrations with "inf" and "60 s" intervals are applied, which increases the signal-to-noise ratio of the continuum image by a factor of 5. The molecular line images were produced by the clean task within CASA with natural weighting. The synthesized beams are 0184 × 0119, 0201 × 0124, 0180 × 0119, and 0202 × 0125, respectively, for SO, SO₂, HCN, and HCO⁺, with position angles (PAs) of 101, 195, 114, and 195. The rms noise levels (σ) for the SO, SO₂, HCN, and HCO⁺ images are 3.5, 4.45, 4.15, and 4.17 mJy beam⁻¹, respectively. The moment maps of the molecular lines were generated by the immoments task using a threshold of 6σ. The continuum image was produced by using line-free channels with uniform weighting. The synthesized beam is 014 × 008, with a PA of 43, and the rms noise level is 0.15 mJy beam⁻¹.

The ALMA observation (Figures 1, 2, and 3) of IRAS 04239+2436 reveals two compact continuum sources corresponding to Sources A and B, without any continuum emission associated with the circumbinary material or beyond the circumbinary scale (∼100 au). Unlike the continuum image, the molecular line emission shows an extended gas structure up to 400 au; the HCO⁺ emission distributes mainly along the bipolar outflow cavities (see Figure 4), while the SO line emission distributes along several arm structures (Figures 1 and 3). The most prominent features are the large extended (∼400 au) multiple-arm (at least three) structures traced by the SO line emission. In contrast, the HCN and SO₂ emissions are detected only toward Source B (Figure 2), which is known as a secondary in the NIR images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The large SO molecular spiral arms probably connect the envelope directly to the circumstellar disks around each stellar component. The arms can be developed by the gravitational torque induced by the dynamical interactions of the binary or multiple protostars (Offner et al. 2010; Matsumoto et al. 2015b). This dynamical process generates shocks, and molecular gas can be excited by shocks to produce detectable emission along spiral arms. A recent study of the sulfur chemistry induced by the accretion shocks (van Gelder et al. 2021) demonstrated that even a low shock velocity (∼3 km s⁻¹) in dense (n ∼ 10⁷ cm⁻³) gas can greatly enhance the abundances of SO and SO₂ via the chemical process initiated by the desorption of CH₄. In the physical conditions developed by low-velocity shocks, H₂S or OCS can also be sublimated from grain surfaces and subsequently react with H, OH, and O₂ to form SO (Millar 1993; Esplugues et al. 2014). Therefore, the SO emission detected along the large arm structures in IRAS 04239+2436 is likely originated by shocks generated by gravitational interaction between the infalling envelope and the orbital motion of the triple system (see Section 3). In contrast, the large dusty arms are not detected around IRAS 04239+2436, probably because of the low dust column density along the spiral arms at scales of several hundreds of astronomical units in the envelope.

The arm structures confined within dense circumbinary disks toward several binary or multiple systems have been detected in the dust continuum (Alves et al. 2019; Diaz-Rodriguez et al. 2022). Theoretical simulations (Bate & Bonnell 1997; Matsumoto et al. 2019b) of a binary system interacting only with a circumbinary disk have shown that two main arms, one of which is connected to each stellar component, are developed, with the secondary component accreting more mass to result in an equal-mass binary at the end. However, IRAS 04239+2436 is likely a triple protostellar system that presents very clear triple arms. In fact, the spiral arms developed by triples in a turbulent envelope (Matsumoto et al. 2015b) have more branches with diverse curvatures compared to the spiral arms developed inside the circumbinary disk (Matsumoto et al. 2019b). Therefore, in order to explain the observed multiple-arm structures in IRAS 04239+2436, we have reanalyzed a hydrodynamic (HD) simulation of a fragmenting turbulent filamentary cloud that forms a multiple stellar system (Matsumoto et al. 2015b). The adaptive mesh refinement (AMR) technique can zoom into one core where multiple protostars are forming (Figure 5).

Zoom In Zoom Out Reset image size

As an initial condition of the simulation, we consider a filamentary cloud in an equilibrium state, in which the thermal pressure supports the cloud against self-gravity in the computational domain of (1.56 pc)³. When we assume isothermal gas of T = 10 K and an infinite cloud length, the density distribution is given by $\rho {\left(R)={\rho }_{0}{\left[1+(R/{R}_{0}\right)}^{2}\right]}^{-2}$ (Stodólkiewicz 1963; Ostriker 1964), where R denotes the cylindrical radius. The scale height of the filament, R₀ = 0.05 pc, is chosen based on the observation from the Hershel survey of the filamentary cloud (Arzoumanian et al. 2011). The density on the filamentary axis is given by ${\rho }_{0}=2{c}_{s}^{2}/(\pi {{GR}}_{0}^{2})=1.45\,\times \,{10}^{-19}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (the corresponding number density is n₀ = 3.79 × 10⁴cm⁻³), where  the isothermal sound speed is c_s = 0.19 km s⁻¹ for gas of 10 K. We impose turbulence with a power spectrum of P(k) ∝ k⁻⁴ as the initial condition, according to the observed size–linewidth relation of σ(λ) ∝ λ^1/2 (Larson 1981), where k, σ, and λ are a wavenumber, the line width, and the length scale, respectively. The average Mach number of turbulence over the computational domain is unity.

The calculation was performed by the AMR code SFUMATO (Matsumoto 2007), assuming the barotropic equation of state, $P{(\rho )={c}_{s}^{2}\rho [1+(\rho /{\rho }_{\mathrm{cr}})}^{7/5}]$, where the critical density is set at ρ_cr = 10⁻¹³ g cm⁻³ (Masunaga et al. 1998). The self-gravity of the gas is taken into account, and the magnetic fields are ignored for the HD simulations. Using the AMR technique, the grid is refined during the course of the calculation, so that the Jeans length of a dense region is resolved by at least eight cells (Truelove et al. 1997). In order to mimic protostar formation, a sink particle is introduced where the density reaches ρ_sink = 10⁻¹¹ g cm⁻³ and several conditions for particle creation are satisfied (Matsumoto et al. 2015a). Each sink particle accretes the gas around it within the sink radius, 2.45 au, increasing the mass.

The evolution of the cloud is as follows. The turbulence disturbs the filamentary cloud, which fragments into cloud cores, because of the gravitational instability (Figure 5(a)). The cloud cores undergo gravitational collapse (Figures 5(b)–(c)); at the center of this, fragmentation produces four protostars (sink particles) in total; and, eventually, three protostars survive (Figure 5(e)). Figure 5 shows a snapshot at t_p = 2.41 × 10⁴ yr, where t_p is the time that has elapsed since the first protostar formation. At this stage, the three protostars exhibit stable orbits, with a close pair with a separation of 20 au (Figure 5(f)) and an outer protostar with a wide separation of 120–200 au. As shown in Figure 5(d), the gravitational interaction between the protostars and the infalling envelope gas produces long spiral arms on a scale of several hundreds of astronomical units.

In this simulation, all the protostars except for the first one are formed by disk fragmentation. Figure 6 shows the time sequence for the formation of protostars and the evolution of their orbits. The first protostar forms in a filamentary structure caused by the turbulence within the cloud core (Figure 6(a)). The protostar has a circumstellar disk and the disk fragments to form a second protostar (Figure 6(b)). Then the third and fourth protostars form in succession, one of which merges with the other, and eventually three protostars survive (Figures 6(c)–(e)). The triple stars show chaotic orbits in the early stages, but the orbits gradually stabilize (Figures 6(f)–(h)). The stage of the simulation when the overall spiral arm structures are matched well with the observation is t_p = 1.95 × 10⁴ yr.

Zoom In Zoom Out Reset image size

The protostars are assigned an ID number in the order of formation. Protostar 1 has the most massive disk, with the highest mass and accretion rate (Figure 7). In contrast, protostars 0 and 2 have low accretion rates and smaller masses. Protostars 0 and 2 form a close binary. Due to a small separation between them, the disks are truncated by dynamical encounters intermittently. These encounters bring about small disks and low accretion rates and also the misalignment of the disks, as shown in Figure 8 (Bate et al. 2010; Bate 2018), which has been considered as a resulting phenomenon of turbulent fragmentation in the early evolutionary stage of binary formation (Lee et al. 2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The interaction between the triplets and the infalling envelope excites the arms. The velocity structure of the arms is shown in Figure 9. In Figure 9(b), the diffuse arms shown to the upper right of the figure have positive radial velocity (outward velocity), while the prominent arms on the left side have negative radial velocity (infall velocity). The infall velocity along the arms reaches 1–1.5 km s⁻¹, and the arm structure can be described as infalling streamers, which are narrow structures that asymmetrically feed gas from the envelope to the scale of the circumstellar disks (e.g., Pineda et al. 2020; Bianchi et al. 2022; Thieme et al. 2022). The spiral arms also have rotation velocity. The arm extending from the south of the triplets has a high rotation velocity due to the gravity of the protostar (Figure 9(c)). Quantitative velocity analysis along the arm is described in Section 4.2.

Zoom In Zoom Out Reset image size

Figure 10 shows the shock velocity and density associated with the arms. When the arms are accelerated by the gravitational torque, the shock velocity exceeds 3 km s⁻¹ (Figure 10(a)). Because an arm accelerated to ∼1 km s⁻¹ collides with the infalling gas with ∼1 km s⁻¹ in the envelope, the shock velocity may reach 3 km s⁻¹. While the arms propagate into the infalling envelope, the shock velocity decreases to 1.5–2 km s⁻¹. The number density in the arms exceeds 10^7–8 cm⁻³ (Figure 10(b)). Such velocities and densities of the arms can enhance the SO abundance along the arms (van Gelder et al. 2021).

Zoom In Zoom Out Reset image size

4.1. A Triple Protostellar System, IRAS 04239+2436

The mid-IR light curve of IRAS 04239+2436 shows a periodicity of ∼8 yr (Figure 11), which is further strong evidence of unresolved close binary interaction within this system, and thus indicative of triple protostars. The anticipated orbital period of the resolved binary is about 270 yr. The finely collimated strong jet and series of jet knots of Source A in IRAS 04239+2436 invoke the existence of disks fueled from a bigger mass reservoir, i.e., the circumbinary disk or envelope (Jørgensen et al. 2022), despite the small sizes of the circumstellar disks truncated by the very close binary interaction within Source A (Reipurth et al. 2000; Reipurth 2000). In addition, the time intervals between the observed jet knots are consistent with the unresolved close binary orbit (Reipurth et al. 2000).

Zoom In Zoom Out Reset image size

In the NIR images, Source A, which is likely an unresolved very close binary, as mentioned above, is brighter than Source B. However, the comparison with the numerical simulation suggests that Source B is the more massive component (Figures 7 and 12), with a larger disk (Figure 8), which can obscure the NIR emission effectively. Another triple system with spirals, L1448 IRS3B (Tobin et al. 2016), also shows that the dimmer component is more massive; unlike IRAS 04239+2436, the more massive one is the close binary in L1448 IRS3B. The HCN line emission is detected only toward Source B, consistently supporting the existence of a very massive disk around Source B (Figure 2(b)), which is likely actively accreting gas along spiral arm 2 (Figures 3 and 12(a)). According to the simulation (see Figure 8), the column density of the larger disk is higher, resulting in its mass being higher than the other two close disks by a factor of 100. In addition, the accretion rate in the primary is higher than those of the secondary and tertiary by more than a factor of 3. This higher accretion rate to the more massive protostar heats its disk more efficiently. As a result, the HCN line intensity, which is proportional to the column density and temperature, is strong enough to be detected in Source B.

Zoom In Zoom Out Reset image size

This accretion streamer associated with Source B and arm 2 is traced by another shock tracer, the SO₂ emission (Figures 2(c)–(d)), as well as the broad SO line profiles around Source B (see Figure 1(b) and the broadened line profiles in Figure 3). In contrast, the 857 μm continuum fluxes toward Sources A (51.5 ± 0.3 mJy) and B (52.8 ± 0.3 mJy) are comparable. This could be due to the continuum optical depth. In the massive disk around Source B, grains easily grow to millimeter sizes and become optically thick at millimeter wavelengths (Harsono et al. 2018; Lee et al. 2019). As a result, the submillimeter emission may not trace the total disk mass, as shown in the ALMA observation of L1448 IRS3B (Tobin et al. 2016).

4.2. Comparison between Observations and Simulations

While the total protostellar mass and their separation are greater than those observed toward IRAS 04239+2436 (Lee et al. 2016), we emphasize that the overall arm features (Figure 12) and their kinematics (Figure 13) produced by the simulation are very similar to the observed ones. Figure 13 compares the kinematics between observations (solid lines) and simulations (dashed lines) along the spiral arms to show that the simulated dynamical process governs the observed arm structures. In Figure 12, the total gas distribution is presented for the simulation, while the observed arm structures in IRAS 04239+2436 are revealed by the SO emission, which traces the shocked gas. The peak intensities of the continuum in the circumstellar disks are about 40 mJy beam⁻¹, and the rms noise level is about 0.1 mJy beam⁻¹. In order to detect the spiral arms in the dust continuum, the column density of the spiral arms must be greater than 1% of the peak column density of the circumstellar disks. While the column densities of the simulated spiral arms are lower than that, they satisfy well the physical conditions for the SO chemistry (van Gelder et al. 2021), as seen in Figure 10.

Zoom In Zoom Out Reset image size

Figure 13(a) presents the velocity at the peak intensity position (hereafter, the peak velocity) along each spiral arm feature in the SO emission as a function of deprojected radius (the symbols connected by solid colored lines). The source velocity is V_LSR = 6.5 km s⁻¹, adopted from C¹⁸O 1–0 and C¹⁷O 1–0 observations (Fuller & Ladd 2002). The peak intensity positions were measured on the peak intensity map of SO 8₈–7₇, which is deprojected by an inclination angle of 55° and a PA of 140° (Figure 13(b)). Because a circumbinary disk has not been detected, we assume 55° and 140°, based on the inclination (50°–61°) and PA (137°–145°) measured from the continuum sources. We also assume that all the spiral arms are on the same plane. For a deprojected radius from 01 to 27, in steps of 01, an intensity profile along the azimuthal direction was extracted and the peak intensity position was measured. At each deprojected radius, the line-of-sight velocity (V_los) was derived from the SO 8₈–7₇ data cube. In Figure 13(b), the peak intensity positions of each arm are marked by symbols in different colors that represent their velocities as given by a color bar. The center position is the center of mass of Sources A and B, on the assumption of equal mass. As seen in the figure, the peak positions mostly well trace the spiral arm features, but arms 2 and 3 are not clearly distinguishable in the inner region (r ≲ 07 or r ≲ 100 au) where SO emission is the brightest, and the two arms are likely overlapped or colliding. Considering the complex emission around the center and the angular resolution (∼01) of the observation, the peak velocity at r < 70 au may not be reliable.

Figure 13(a) also presents the velocities of the arms derived from our simulation (dashed colored lines) and the Keplerian rotation (black lines) for comparison. Since the total mass and separations among the companions from the simulation are larger than those of IRAS 04239+2436 by a factor of 2, the deprojected radii of the simulation are scaled by a half in the figure. As the overall arm features of IRAS 04239+2436 and the simulation are similar (Figure 12), their velocity profiles also show similar trends (Figure 13(a)). The velocity of the arms is not explained by a single Keplerian rotation curve, although arms 1 and 2, beyond 100 au, roughly trace the Keplerian rotation caused by the mass between 0.6 M_⊙ and 0.15 M_⊙. In addition, the observed structure of arm 1 beyond 15 is not fitted by a polynomial function simultaneously with the inner part of the arm as presented in Figure 13(b).

Note that Figure 13(a) displays the line-of-sight velocity, which is a composite of the infall velocity and the rotation velocity. In contrast, Figure 14(a) shows the infall velocity and rotation velocity separately along the arms, based on the simulation with 3D velocity components. The location of the arms in the simulation is depicted in Figure 14(b). Arm 1 (represented by the red lines in Figure 14(a)) has a considerable infall velocity, which is comparable to, but less than, the rotation velocity in the range of 150 au ≲ R ≲ 300 au. Arms 2 and 3 (blue and green lines) exhibit considerably larger infall velocity than rotation velocity in the outer part where R ≳ 100 au. Therefore, arms 1, 2, and 3 are classified as infalling streamers, but their infall velocities are less than the freefall velocity (black solid line). In the inner region where R ≲ 100 au, all the arms exhibit rotation velocities that exceed the infall velocities and the Keplerian velocity (black dashed line) due to the gravitational torque from the orbiting protostars.

Zoom In Zoom Out Reset image size

4.3. Multiple Formation and Magnetic Field Strength

A magnetic field must be weak to form a multiple stellar system via the fragmentation of a cloud core or a disk (Machida et al. 2005). To test the effect of the magnetic field on the formation of multiples and the following development of large spiral arms, we also ran magnetohydrodynamic (MHD) simulations, varying the initial magnetic field strength, and present the results in Figure 15.

Zoom In Zoom Out Reset image size

The initial condition and numerical models are the same as those of the HD model, but the magnetic fields are included. Uniform magnetic fields perpendicular to the filamentary cloud are imposed as the initial condition, because recent observations have suggested that dense filaments tend to have perpendicular magnetic fields (Pattle et al. 2022). Three models with the initial magnetic field strengths of 1, 5, and 10 μG are considered. The magnetic field strength is often described in terms of the mass-to-flux ratio, and the critical magnetic field strength (Nakano & Nakamura 1978; Tomisaka et al. 1988) is estimated as B_cr = 2π G^1/2Σ = 20 μG for our models. The initial magnetic field strengths of 1, 5, and 10 μG are therefore 0.05, 0.25, and 0.5B_cr, respectively (corresponding to the nondimensional mass-to-flux ratios of μ = 20, 4, and 2). Thus, only magnetically supercritical clouds are considered here. For an MHD scheme, we adopted Boris-HLLD (Matsumoto et al. 2019a), which allows us to follow a long-term evolution, even in cases of strong magnetic fields. The ohmic dissipation is considered according to previous simulations (Matsumoto 2011; Matsumoto et al. 2017). The MHD scheme has second-order accuracy in space to maintain numerical stability, while the HD scheme has third-order accuracy.

The MHD models exhibit similar evolution to the HD model; the filamentary cloud fragments into cloud cores, which undergo gravitational collapse. A multiple- or single-star system forms at the center of the cloud core, depending on the initial magnetic field strength. Models with weaker magnetic fields tend to form multiple systems with more companions (Figure 15) and to have more prominent and more extended arms (Figure 16). The models with 0, 1, 5, and 10 μG form the systems consisting of three (triple), four (quadruple), two (binary), and one (single) stars, respectively.

Zoom In Zoom Out Reset image size

According to our MHD simulations, the natal cloud of IRAS 04239+2436 must have a very weak magnetic field, with a nondimensional mass-to-flux ratio larger than μ = 2–4, to form a multiple system. Recent observations of the submillimeter dust emission polarization also show that the strong magnetic field in large filamentary structures is weakened down to the magnetically supercritical condition at small core scales (Doi et al. 2021). The present observation provides more strong observational evidence of the reduced effect of the magnetic field in a star-forming dense core (Eswaraiah et al. 2021).

4.4. Planet Formation in Multiple Protostellar Systems

Although most stars form in multiple-star systems, their formation process is still controversial. In contrast to the single-star formation process, dynamical interactions take place not only among the protostars, but also between the central compact objects and the infalling gaseous envelope. Because the accretion history can be lost in a few orbital timescales, it is important to capture the signatures of these dynamical interactions by zooming into such a system at the earliest possible formation epoch in order to better understand the formation process of multiple-star systems. According to our high-resolution and high-sensitivity ALMA observation, IRAS 04239+2436, which has been known as a protostellar binary system, shows clear triple spiral arms surrounding a young potentially triple protostellar system in the SO line emission. The spiral arms traced by shocked SO gas act as accretion flows over 400 au in length from the large-scale envelope to the 50 au scale stellar system. The imaging results are compared with the numerical MHD simulations of analogous systems. While both magnetic fields and turbulence are important to shaping the natal molecular cloud, our numerical simulations suggest that the multiple asymmetric arms may be characteristic of multiple protostars forming in regions of weak magnetic fields.

ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.01397.S. J.E.L. is supported by a National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT; grant No. 2021R1A2C1011718). D.H. is supported by the Centre for Informatics and Computation in Astronomy (CICA) and grant No. 110J0353I9 from the Ministry of Education of Taiwan. D.H. acknowledges support from the Ministry of Science of Technology of Taiwan through grant No. 111B3005191. Numerical computations were carried out in part on XC50 (ATERUI II) at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. T.M. is supported by JSPS KAKENHI grant Nos. JP23K03464, JP18H05437, and JP17K05394. This publication also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration. K.T. was supported by JSPS KAKENHI grant No. JP20H05645.

Facility: ALMA - Atacama Large Millimeter Array.

Software: CASA (McMullin et al. 2007), SFUMATO (Matsumoto 2007).