The First Bird's-eye View of a Gravitationally Unstable Accretion Disk in High-mass Star Formation

Figure 1(a) shows the dust continuum image of the accretion disk that consists of the compact central emission and resolved structure. The latter was highlighted in Figure 1(b), where the former compact emission was subtracted by elliptical Gaussian fitting (Appendix A). The compact emission has an averaged diameter of ∼160 au (Table 1). The peak T_B is 480 K. The resolved structure is almost "ring-like" with the outer radius of 250 au. There are an arc-like elongated feature in the east and several clumpy features in the west. The peak T_B in the east side (up to 126 K) is 2–3 times brighter compared to other, less bright parts. This asymmetry would reflect the contrast of the surface density rather than the asymmetric temperature profile, as the disk temperature should be basically dominated by stellar radiation in this small scale. It should be noted that, if the disk actually has a non-axisymmetric structure such as spiral arms, the dynamical heating effect can cause a deviation from the axisymmetric temperature profile. Although few limited studies observed the innermost accretion system in high-mass star formation (Cesaroni et al. 2014; Beuther et al. 2017; Hirota et al. 2017), this is the first completely resolved view of such a compact non-axisymmetric structure inside a disk.

Zoom In Zoom Out Reset image size

The compact emission is consistent with the continuum emission at 45 GHz previously reported in Motogi et al. (2017) as shown in Figure 1(b) and Table 1, indicating that we could successfully subtract the hot and compact emission in the innermost region. The spectral index is 2.5 ± 0.2 between 45 and 150 GHz, considering the flux uncertainty of 10%. The optical depth, averaged surface density (N_H2), and total mass of the compact emission are ${2.0}_{-0.9}^{+1.4}$, ${8.2}_{-3.8}^{+5.7}\times {10}^{25}\,{\mathrm{cm}}^{-2}$ and ${0.8}_{-0.4}^{+0.5}\,{M}_{\odot }$, respectively (Appendix A). On the other hand, the optical depth in the resolved structure was determined as ${0.68}_{-0.30}^{+0.36}$ at the eastern peak and 0.31 ± 0.12 on average, which corresponds to the surface densities of ${2.8}_{-1.2}^{+1.5}\times {10}^{25}\,{\mathrm{cm}}^{-2}$ and 1.3 ± 0.5 × 10²⁵ cm⁻², respectively. The total mass of the dusty accretion system is ${1.9}_{-0.8}^{+1.0}\,{M}_{\odot }$, including the compact component. Here we adopted the dust mass opacity of 0.59 cm² g⁻¹ at 150 GHz and the gas-to-dust ratio of 100 (Appendix A).

Figure 1(c) shows the entire continuum emission in a further large-scale envelope (∼700 au in radius). The integrated flux density is 316 mJy (>3σ level), including 140 mJy from the compact emission within 250 au. The mass outside of 250 au is ∼5 M_☉, which is estimated by the graybody emission adopting the same dust parameters and the averaged temperature of ∼100 K that is deduced from line emissions (see the next subsection). The total mass of the entire system within 700 au radius is 7 M_⊙. The uncertainty of these mass estimations originates primarily from the uncertainty of the dust parameter. This will be discussed later.

We have detected seven CH₃OH lines (Figure 2(a) and Table 2). Figure 2(b) presents the integrated intensity map of the molecular envelope that has a comparable scale with that of the entire dust emission. Here, we stacked three CH₃OH lines, with the upper-state rotational energy higher than 70 K (Appendix B). It should be noted that all of the detected lines are thought to be not maser but thermal emission, because of the extended structure and moderate brightness. The envelope size is 700 au in radius and completely encloses the compact dusty structure. The CH₃OH emission has significantly decreased toward the center, which is consistent with the optically thick dust emission. Almost all of the CH₃OH emission is from beyond 250 au in radius. Figure 2(c) shows the position–velocity (PV) diagram of the CH₃OH emission along the north–south direction (Appendix B). The entire velocity structure is roughly symmetric against the systemic velocity of the natal cloud (∼−4.5 km s⁻¹).

Zoom In Zoom Out Reset image size

We detected both the rotational spin-up motion and infalling motion at the same time (e.g., Sakai et al. 2014; Beuther et al. 2017). These kinematic features could be basically reproduced by a simple model of the infalling rotating envelope, as in the case of low-mass protostars (Sakai et al. 2014; Oya et al. 2016). Figure 2(d) shows the schematic view of our model, which is described in Appendix B. The envelope is rotating counterclockwise around the rotational axis that is inclined 8° from the LOS (Motogi et al. 2016, 2017). The envelope becomes thicker outward with a constant flaring angle. We adopted the specific angular momentum of 2.4 × 10²¹ cm² s⁻¹ in the model. This corresponds to the centrifugal radius of ∼250 au, which is consistent with the outer radius of the resolved structure in Figure 1(b) assuming the stellar mass of 10 M_☉.

We calculated the LOS velocities on the midplane and the near/far side of the envelope surface (Figure 2(d)). We have found that the envelope should have a relatively wide flaring angle (∼35°), in order to reproduce the rotational spin-up motion and the centrifugal brake against the infall acceleration. For example, if the envelope is geometrically thin with a much smaller flaring angle, the velocity field coincides with that in the midplane. On the other hand, the specific angular momentum that is significantly smaller than the assumed value could not reproduce the feature of the centrifugal braking. Conversely, the rotational spin-up feature at the angular offset of ±02 could not be reproduced by the significantly larger angular momentum. Therefore, we suggest that the inner dust emission traces a centrifugal disk and the rotating accretion flow reaches Keplerian rotation inside 250 au. It should be noted that further fine-tuned model fitting is beyond the scope of this study, because we clearly require higher spectral resolution and better sensitivity to detect any faint high-velocity emission.

We also detected four ¹³CH₃OH lines (Figure 2(a) and Table 2). The spatial structures of those ¹³CH₃OH lines are almost identical to that of the main isotopes. Surprisingly, the observed line ratio (CH₃OH/¹³CH₃OH) is very small (∼1.3), suggesting that these line emissions are optically thick and almost thermalized. The averaged T_B of these optically thick CH₃OH lines is ∼100 K. We used this value as the dust temperature in the outer envelope.

If we adopt the isotopic abundance ratio ¹²C/¹³C ∼ 60, the averaged optical depth of ∼90 is expected for the main isotopes. By using the non-local thermodynamic equilibrium radiative transfer code (RADEX: van der Tak et al. 2007), we estimated the lower-limit density (n_H2) and column density (N_CH3OH) as 3 × 10⁹ cm⁻³ and 10¹⁹ cm⁻², respectively. Here, we adopted a kinematic temperature of 100 K and velocity dispersion of 5 km s⁻¹, which is the observed line width (FWHM). The lower-limit N_H2 is ∼10²⁵ cm⁻², assuming the relative CH₃OH abundance of 10⁻⁶, which is the highest abundance observed in hot cores (e.g., Bisschop et al. 2007). We finally obtained a total envelope mass of 7 M_⊙ using the averaged radius of 700 au. This is comparable with that estimated from the dust emission (∼5 M_⊙).

We note that actual dust temperature might be colder than the gas temperature estimated by line emission, because of different spatial locations in the disk/envelope as suggested in recent observational studies (e.g., Beuther et al. 2017). The envelope mass estimated by dust emission could be larger in this case.

The kinematic model indicates that an infalling timescale of the molecular envelope is ∼10³ yr. This was estimated by numerically solving the equation of motion along the envelope surface for the mass particle at the outer edge. Here we only considered the stellar gravity and centrifugal force in the equation. On the other hand, the averaged envelope density of 3.5 × 10⁻¹⁵ g cm⁻³, where 6 M_⊙ is divided by the volume of the modeled envelope, corresponds the freefall time of ∼1.1 × 10³ yr. These are comparable, but the numerical value is slightly shorter than the freefall time because of the stellar gravity. The averaged infall rate inside the envelope is ∼6 × 10⁻³ M_⊙ yr⁻¹. We expect an effective accretion rate onto the stellar surface as 3 × 10⁻³ M_⊙ yr⁻¹ or higher, because 50% of accreting mass could be converted to the outflow in case of high-mass star formation (e.g., Matsushita et al. 2017).

This infall rate is clearly larger than that estimated from the Jeans mass (i.e., $\sim {c}_{s}^{3}/G$). It becomes ∼10⁻⁴ M_⊙ yr⁻¹, even adopting the current envelope temperature of 100 K. This discrepancy is usual in high-mass star formation, where we must consider the so-called effective Jeans mass, including a non-thermal velocity component such as turbulent or Alfvén velocities (e.g., Hosokawa & Omukai 2009; Inoue & Fukui 2013). We can replace c_s with c_eff, which is the root sum square of the velocity components in this case. The observed infall rate corresponds c_eff of ∼3 km s⁻¹, suggesting a highly turbulent and/or strongly magnetized condition of the initial core. Such a condition may be achieved by dynamical compression via cloud–cloud collision (e.g., Inoue & Fukui 2013), rather than simple gravitational condensation.

The age of G353 could be only 3 × 10³ yr, considering the infall rate of 3 × 10⁻³ M_⊙ yr⁻¹. We note that the accretion rate onto the stellar surface from the disk can be highly variable in the case of the gravitationally unstable disk; however, the time-averaged accretion rate should be consistent with the infall rate onto the disk from the envelope (e.g., Meyer et al. 2018) after subtracting the outflow rate. We suggest that G353 is one of the youngest (<10⁴ yr) high-mass objects currently known, even considering an overall error of factor 3. The final stellar mass could be 13–17 M_⊙ depending on the outflow efficiency. It should be noted that this is still the lower-limit mass, as any envelope structure larger than 05 should be resolved out in our data. Further observations on a larger scale are required to confirm the final stellar mass and its evolutionary stage.

The disk mass within 250 au reaches almost 20% of the stellar mass; therefore, the disk could be self-gravitating (Forgan et al. 2016). Figure 3(a) shows the spatial distribution of Toomre's Q parameter (Toomre 1964; see Appendix C), which is a well-known measure of the gravitational stability (Q < 1: the disk immediately collapses; 1 < Q < 2: the disk is unstable against non-axisymmetric perturbations). The central region shows Q > 2, suggesting that the disk is stabilized by the stellar radiation and stronger gravity of the central star. On the other hand, a significant fraction of the outer disk shows that Q < 2 and the minimum Q reaches 0.9 around the eastern local peak. This fact indicates that the disk is, at least, unstable against non-axisymmetric perturbations such as a spiral arm (e.g., Tomida et al. 2017). Furthermore, the averaged Q of ∼1.9 outside 80 au in radius suggests that the entire disk is moderately unstable, avoiding immediate collapse during the freefall time (∼a few hundred yr). Such a short timescale is too short to be observed.

Zoom In Zoom Out Reset image size

Although these criteria are valid for low-mass star formation, it may not be the case for high-mass star formation. Recent high-resolution numerical studies in Meyer et al. (2018) pointed out refined criteria (Q < 1 for spiral arms, Q < 0.6 for fragmentation). They also found that a disk-to-star mass ratio is 0.5 or higher in the case of the unstable disk. These new criteria suggest that our observed Q could be a little too high for the disks to be gravitationally unstable. The major uncertainty of Q in this study originates from the uncertainty of the dust opacity index β, which is used for estimating the surface density Σ. The uncertainty from errors in the sound velocity c_s and angular velocity Ω is limited, because these are proportional to the square root of a temperature or stellar mass. We note that β could vary within 1–2, as discussed in Appendix A. Although we adopted β of unity for conservative mass estimation; if β = 2 is adopted, the dust mass opacity κ_ν becomes 40% smaller. The total disk mass reaches 7 M_⊙ within 250 au (5 M_⊙ for the compact source within 100 au) in this case. We also obtain a 13% lower dust temperature, because higher optical depth is expected at the center (see Appendix A). This results in the averaged column density of 5.3 × 10²⁶ cm⁻² for the compact source and 2.4 × 10²⁵ cm⁻² for the resolved disk. Recalculated Q values are 0.9 (average), 0.4 (minimum), and 2.0 (maximum) in the resolved disk. This is more consistent with the results of Meyer et al. (2018).

The spatial scale of the simulated disk in Meyer et al. (2018) is quite similar to that of our observed disk, allowing us to directly compare disk structures. For example, we estimated the gas density of the disk, using the observed column density and disk scale height (H ∼ c_s/Ω). We obtained H ∼ 13 au, at the radius of 100 au (see Appendix C for c_s and Ω estimation). The gas density at 100 au is estimated as 1.7 × 10⁻¹² g cm⁻³ for β = 1 and 1.2 × 10⁻¹¹ g cm⁻³ for β = 2. The latter value is comparable to that reported in Meyer et al. (2018). These comparisons could imply that β = 2 is a more adequate choice for G353; this will be examined by ongoing follow-up multi-band ALMA observations.

Although β uncertainty remains, the disk in G353 is the most self-gravitating disk currently known. We propose that this is the first detection of a gravitationally unstable accretion disk in high-mass star formation. The obtained Q values are significantly smaller than those estimated in previous studies toward high-mass star-forming regions (Beuther et al. 2017; Ahmadi et al. 2018). This fact could be also explained by the young age of G353; i.e., the host object is still relatively less massive compared to other disk candidates, and the disk itself consists of gas with a smaller angular momentum at the inner region of the natal core. It is evident that radiative heating from the 10 M_⊙ accreting star is insufficient to stabilize the young massive disk with the radius of 250 au. The highly time-dependent outflow activity (Motogi et al. 2011, 2016) may be related to the unstable nature of the accretion flow.

The detected asymmetric structure can be also naturally explained by self-gravitating effects. The observed surface density contrast (up to 4–5; Figure 3(b)) is easily produced by a spiral arm. If this is the case, the gravitational torque can redistribute angular momentum in the disk. Such an effect can produce further infalling streams within the inner region (<100 au) found in the previous study (Motogi et al. 2017). The specific angular momentum of 2.4 × 10²¹ cm² s⁻¹ in the envelope is four times larger than that estimated for the inner streams. This indicates that 70% of the angular momentum is removed in the disk, probably via the gravitational torque. A rapid mass accretion to the central star promoted by non-axisymmetric structures decreases the disk surface density; therefore, Q should be adjusted to Q ∼ 1–2 (e.g., Tomida et al. 2017) as seen in Figure 1(d).

Another possible origin of such asymmetric structures is a tidal perturbation via a close-encounter event by a nearby cluster member (e.g., Thies et al. 2010). However, this is unlikely to be the case because no bright cluster member was detected in the ALMA FOV. Our continuum sensitivity corresponds to a 1σ dust mass sensitivity of 8 × 10⁻³ M_⊙ beam⁻¹, adopting a dust temperature of 100 K. It is clear that G353 is the first high-mass cluster member in the natal cluster-forming clump, although there may be a very small dust condensation that could form a brown dwarf.

The clumpy local peaks seen in the western region could also trace the fragmentation process in the disk. Indeed, the outer disk radius of ∼250 au and the minimum Q parameter of ∼1 are remarkably similar to that in the fragmented disk with the triplet low-mass protostar system in L1448 IRS3B (Tobin et al. 2016). These clumps may be the seed of a lower-mass binary companion, or accrete on the central object causing the accretion burst event (e.g., Meyer et al. 2017, 2018, 2019). The latter event was recently detected in several objects (Caratti o Garatti et al. 2017; Hunter et al. 2018). The masses of two large clumps are 0.04 and 0.09 M_⊙ (Appendix D). These masses are for the β = 1 case, and they are twice as large for the β = 2 case.

For comparison, we estimated the Jeans scale (λ_J) and Jeans mass (M_J), in the disk. For the β = 1 condition, the gas density is 1.6 × 10⁻¹³ g cm⁻³ and the temperature is ∼230 K at 250 au radius. The former was estimated from the averaged surface density of 1.0 × 10²⁵ cm⁻² in the eastern region, divided by the scale height of ∼43 au. We finally obtained λ_J of ∼104 au (006) and M_J of ∼0.07 M_⊙. This is consistent with the observationally estimated masses of the clumps, and is still comparable even for the β = 2 condition, where λ_J and M_J decrease by 40%. Although the estimated clump masses are 10–20 times larger than the accreting mass in the burst event in S255IR (Caratti o Garatti et al. 2017), it is consistent with the clump mass range formed by the disk fragmentation events in the numerical studies in Meyer et al. (2019).

Further high-resolution (<005) ALMA observations are required in order to directly investigate the innermost region. In this case, it is possible that the dust emission could be completely optically thick in a higher frequency than ALMA band 4. Lower-frequency interferometers such as the Square Kilometer Array or Next Generation Very Large Array may be able to resolve the innermost accretion process.