A Hot and Massive Accretion Disk around the High-Mass Protostar IRAS 20126+4104

What role accretion disks play in the formation of high-mass stars ($M\gtrsim 8M_{\odot }$) remains a long-standing question. Circumstellar disks form naturally in the centers of rotating inflows and are a key element in the standard paradigm of the formation of Sun-like stars, providing for the growth of planetary systems. However, it remains debatable whether high-mass stars form in a similar fashion. The accretion rates in massive star formation may be high enough to induce gravitational instabilities and put disk-mediated accretion in doubt. Previous observations reporting disk-like accretion flows around high-mass protostars have not yet assessed the stability of the candidate disks (Chini et al. 2004; Patel et al. 2005; Jiménez-Serra et al. 2007; Cesaroni et al. 2014; Johnston et al. 2015). On theoretical grounds, a very high accretion rate ($\gtrsim {10}^{-4}{M}_{\odot }y{r}^{-1}$) is required to form stars more massive than $8M_{\odot }$. Stars this massive undergo rapid enough Kelvin-Helmholtz contraction that they begin hydrogen burning while still accreting (Palla & Stahler 1993). A continuous resupply of fresh hydrogen is required to allow a growing massive protostar to reach the mass of a B or O star before exhausting its hydrogen fuel and leaving the main sequence (Keto 2003). The gas densities required for such rapid accretion could induce gravitational instabilities in the accretion flow that would make the disk prone to fragmentation, perhaps producing companion objects (Kratter & Matzner 2006). On the other hand, the stability of a disk is determined not only by the gas density, but also the differential shear and the gas temperature which are much affected by the disk environment (Durisen et al. 2007). High stellar luminosity ($L>{10}^4{L}_{\odot }$) and shock heating may keep the disk warm enough to allow accretion to proceed steadily (Pickett et al. 2000b; Kratter et al. 2010). Massive disks are also subject to dynamical heating from spiral shocks, and vertical shear at the interface with the envelope (Pickett et al. 2000a; Harsono et al. 2011). The stability of the disk may be assessed by the Toomre-$Q$ parameter,

which compares the stabilizing effects of the temperature (sound speed $c_s$) and shear (angular velocity $\Omega$) against the clumping tendency toward instability induced by the surface density ($\Sigma$). Values of $Q>1$ imply a stable disk. Up until now, there has been insufficient observational evidence to answer the question. We present new molecular line observations of the massive protostar IRAS 20126+4104 that directly address the question of the gravitational stability of massive protostellar disks.

IRAS 20126+4104 (hereafter I20126) is a nearby ($1.64kpc$; Moscadelli et al. 2011), luminous ($\sim 1.3\times {10}^4{L}_{\odot }$; Johnston et al. 2011) high-mass protostar with a well collimated bipolar molecular outflow (Zhang et al. 1999; Moscadelli et al. 2005; Su et al. 2007; Hofner et al. 2007) and an accretion disk (Zhang et al. 1998; Cesaroni et al. 2005, 2014; Xu et al. 2012). Previous radio and infrared observations suggest a disk mass of several $M_{\odot }$ and a protostellar mass of $7$ to $12M_{\odot }$ (Sridharan et al. 2005; Cesaroni et al. 2005; Keto & Zhang 2010; Johnston et al. 2011). With such a high ratio of the disk to stellar mass, the self-gravity of the disk is significant in the accretion dynamics, and the stability of the disk is questionable. A previous study compared observed infrared emission against that predicted by a model of a disk-mediated accretion flow (Johnston et al. 2011) and found the Toomre-$Q$ high enough to allow self-consistency with their original assumption of a disk. Our new observations directly measure the temperature, angular velocity, and surface density to enable a direct determination of the Toomre-$Q$ and the disk stability.

Our observations of the $C{H}_{3}CN$ line emissions around I20126 with the Submillimeter Array⁴⁴4The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica. (SMA; Ho et al. 2004) achieved the highest possible angular resolution ($0\text{farcs}37$, equivalent to $\sim 600AU$) sufficient to spatially resolve the accretion flow. The $C{H}_{3}CN$ molecule is ideal for the identification of the hot, dense gas in the flow near the star. It requires high gas densities, $\gtrsim {10}^{4}c{m}^{-3}$, for collisional excitation, and high temperatures, $\gtrsim 100K$, to produce detectable emission (Araya et al. 2005; Chen et al. 2006). The Doppler shifting of the molecular lines measures the rotational velocities and shear within the accretion disk. The emission includes multiple line transitions with different excitation energies allowing a measurement of the gas density from the line brightness and a measurement of the gas temperature from the brightness ratios.

The present study of the high-mass protostar I20126 finds its massive disk hot and stable against gravitational fragmentation even as the accretion proceeds at a high rate. This paper is organized as fellows. Sect. 2 describes details of the observations and data reduction. We present the basic observational results in Sect. 3 and explain details of radiative transfer models for both continuum and line emissions in Sect. 4. We discuss the implication of our results in Sect. 5 and conclude with a brief summary in Sect. 6.

The $870\mu m$ dust emission (Fig. 1) shows a rather symmetric distribution with the peak position at $(\alpha ,\delta )\left(J2000\right)=$ (20:14:26.029,+41:13:32.57), which is assumed to be the position of the protostar. Such symmetric morphology already suggests a disk at a moderate inclination, but the character of a rotating disk is revealed through spectral imaging. The observed spectrum covers the $K=0,1,\dots ,9$ components of the $C{H}_{3}CN$ $J=19-18$ transition with upper state energy in the range of $E_{up}=168-745K$ and the $C{H}_{3}OH$ (${14}_{1,13}-{14}_{0,14}$) emission with $E_{up}=260K$ (Fig. 2). The $C{H}_{3}CN$ lines are collisionally excited and thus sensitive to gas temperature (Araya et al. 2005; Chen et al. 2006). We can compare the Doppler shifting of spectral lines of different excitation energies to identify the spin-up of the accretion flow and the temperature gradient within disk created by the hot protostar. The excitation is indicated by two quantum numbers, $J$ and $K$, which specify the total angular momentum of the molecule and its projection along the principal axis. As examples, we showed the $K=2,3,6,9$ components of the $C{H}_{3}CN$ $J=19-18$ transition (Fig. 3a–3d), whose upper states have increasing energies of $E_{up}=196,232,425,745K$ and are expected to trace emission progressively close to the hot protostar. In addition, the $C{H}_{3}OH$ (${14}_{1,13}-{14}_{0,14}$) line shows a similar and slightly more extended emission, likely to trace part of the envelope (Fig. 3e).

The spectral line data are three-dimensional (3-D) with axes of position, position, and velocity and can be displayed in several ways. Plotting the average spectral line velocity (first moment) at each position shows a consistent gradient characteristic of rotation for all five observed lines (Fig. 3a–3e). We can select a single plane in the data cube oriented along this gradient, and plot the line intensity as a function of position and velocity (P-V diagram). These plots (Fig. 4a–4e) reveal progressively steeper velocity gradients with increasing $E_{up}$, indicating an increase in the rotational velocity of the higher temperature gas nearer the star. Such spin-up motions are characteristics of an accretion flow that at least partially conserves angular momentum, the limit of which is a rotationally supported or Keplerian disk.

The velocity gradients in the spectral line data also show the inward or radial flow of accretion in the more circular and bluer pattern of the average velocity of the $K=2$ line (Fig. 3a) compared with $K=9$ line (Fig. 3d). This effect arises from self-absorption of the spectral line emission in a radial flow. Along a line of sight through the center, a radial flow splits the line emission into red and blue shifted components from the near and far sides of the core. If the line were optically thin, we would see a symmetrically split profile. However, absorption by colder gas in the outer part of the core selectively absorbs the red-shifted emission, which is closer to its own velocity, while the emission from the far side, which is blue-shifted to a dissimilar velocity, passes through. If the flow were purely radial, this would produce a bullseye pattern of velocities, bluest in the center. In Fig. 3a, we see a combination of this effect along with the rotation. In contrast, a purely rotational flow will show a simpler one-dimensional gradient across the image as shown in Fig. 3d. This indicates that the high temperature gas has spun-up so much that the rotational velocities dominate the average. In I20126, the spectral line velocities and brightnesses are consistent with a more radial flow in a cooler accreting envelope that spins up and flattens to a hot rotating disk as the accretion flow approaches the star.

The velocity patterns provide a qualitative overview of the accretion flow feeding I20126. We can extract precise measurements of the disk and envelope temperatures, densities, and velocities by comparing the observations with model accretion flows. Following earlier studies (Keto & Zhang 2010; Johnston et al. 2011) that modeled lower angular resolution observations of other spectral lines, e.g. $N{H}_3$, and the infrared continuum observations of I20126, we constructed a 3-D analytical model describing a thin accretion disk (Pringle 1981) in Keplerian motion enveloped within the centrifugal radius of an angular-momentum-conserving accretion flow (Keto & Zhang 2010; Ulrich 1976). We include the stellar irradiation to heat the flared disk (Kenyon & Hartmann 1987) consistent with the presence of outflow cavities (Qiu et al. 2008; Moscadelli et al. 2011). A bolometric luminosity of $L_{bol}=1.3\times {10}^4{L}_{\odot }$ is assumed (Johnston et al. 2011). We also adopt the systemic velocity of $-3.5km{s}^{-1}$ (Cesaroni et al. 1999).

Under conditions of local thermodynamic equilibrium (LTE), we solved the radiative transfer equation for the intensity, $I_{\nu }$, of the continuum and spectral lines simultaneously to construct synthetic continuum images and 3-D spectral image cube for comparison with the observations. The radiative transfer is evaluated with the source function, i.e. the Planck function, $B_{\nu }\left(T\right)$, modulated by the linear sum of opacities for the continuum and spectral lines. At the center frequency of each individual channel, we solve for $I_{\nu }$ with

where ${\alpha }_{\nu }$ is the total absorption coefficient given by

For the continuum emission, we have ${\alpha }_{\nu }^{cont}=\left(\rho /100\right){\kappa }_{\nu }$, where $\rho$ is the mass density, and a gas-to-dust mass ratio of $100$ is assumed. We use the dust opacity law ${\kappa }_{\nu }=10(\lambda /250\mu m{)}^{-\beta }c{m}^2{g}^{-1}$ (Hildebrand 1983) with $\beta =1.5$ for the disk (Sec. 2.2) and $\beta =1.8$ for the envelope (Johnston et al. 2011). Regarding spectral lines, the line blending is significant for a few molecular lines so we have

where $N_{line}$ is the total number of lines included in the model. The absorption coefficient of the $i$-th line, ${\alpha }_{\nu }^{\left(i\right)}$, is given by

where $n$ is the gas density, $X_{mol}$ the abundance of the molecule, $g_{up}^{\left(i\right)}$ the upper state degeneracy, $E_{up}^{\left(i\right)}$ the upper state energy, $Q_{mol}\left(T\right)$ the partition function of the molecule, $A_{ul}^{\left(i\right)}$ the spontaneous emission rate of the line, ${\nu }_0^{\left(i\right)}$ the rest frequency of the line, and ${\Phi }_{\nu }$ the model line profile (see Appendix A).

The integration along line of sight is performed using the Runge-Kutta method in steps of $0.01-300AU$. For comparison with the observed spectral line data, synthetic spectral line images are generated through observation simulation, including visibility sampling and image making, with MIRIAD. The Levenberg-Marquardt method was used to optimize models by minimizing ${\chi }^2$ value computed with both the continuum and spectral line data. Only channels with emission stronger than $3\sigma$ in the central $2\text{arcsec}$ region (dark gray histogram in Fig. 2) are considered for the model fitting and the ${\chi }^2$ computation. The model spectral imaging includes 19 molecular lines (Table 1): the $K=0,1,\dots ,9$ components of the $C{H}_{3}CN$ $J=19-18$ transition, the $K=0,1,\dots ,7$ components of the $C{H}_3^{13}CN$ $J=19-18$ transition, and $C{H}_{3}OH$ (${14}_{1,13}-{14}_{0,14}$). The ${}^{12}C/{}^{13}C$ ratio is assumed to be 70 for a galactocentric distance of $8.3kpc$ (Wilson & Rood 1994). The model contains nine adjustable parameters as listed in Table 2. The uncertainty of the parameter $a_i$ is estimated by $\Delta a_i=\pm \sqrt{C_{ii}}\sqrt{N_{beam}}\sqrt{\Delta {\chi }^2}$, where $C_{ii}$ is the $i$-th diagonal term of the covariance matrix, $N_{beam}$ is the number of pixels in the synthesized beam, and $\Delta {\chi }^2\equiv {\chi }^2-{\chi }_{min}^2=10.43$ gives the 68.3% confidence level for degrees of freedom of nine.

Analogous to the standard model of star formation with disk-mediated accretion, our model consists of a thin disk (Pringle 1981) of mass $M_d$ in Keplerian motion around a stellar mass, $M_{\ast }$, residing within the centrifugal radius, $R_c$, of an accretion flow with constant specific angular momentum, $\Gamma$, in an infalling envelope (Ulrich 1976). To account for a fairly large disk mass, the centrifugal radius is computed with

Density-weighted means are calculated for velocity and temperature in positions where both disk and envelope are present. Due to the prominent outflow cavities (Qiu et al. 2008), stellar irradiation is also included for anticipated heating due to the flared disk geometry (Kenyon & Hartmann 1987). Considering various heating processes, e.g. accretion shocks, that may occur to further raise the disk temperature, we introduce a scaling factor, $B_T$, to take these effects into account. Due to a fairly large disk mass, corrections for the enclosed disk mass as a function of radius are applied to rotation velocities and viscous heating produced by differential shear in the disk.

The density distribution of the flared disk is described as

where $H\left(R\right)=H_0(R/R_{\ast }{)}^{1.25}$ is the disk scale-height with $H_0=0.01R_{\ast }$, and ${\rho }_{d0}$ is related to the disk mass, $M_d$, by ${\int }_{R_{in}}^{R_c}{\int }_{-\infty }^{\infty }{\rho }_d(R,z)dz2\pi RdR=M_d$, where $R_{in}$ is the inner radius of the disk. The surface density is defined by $\Sigma \left(R\right)\equiv {\int }_{-\infty }^{\infty }{\rho }_{d}dz\propto R^{-1}$, and the enclosed disk mass is given by $M_d\left(R\right)={\int }_{R_{in}}^R\Sigma \left(R^{\prime }\right)2\pi R^{\prime }d{R}^{\prime }$. The Keplerian rotation speed in the disk is hence

Given a disk mass accretion rate, ${\dot{M}}_d$, and the disk surface density, $\Sigma \left(R\right)$, an inward accretion velocity can be computed from

which is roughly constant through the disk.

We describe the central massive protostar as a zero-age main-sequence star (Schaller et al. 1992) with additional surface heating by gas accreted from the inner edge of the disk, releasing all its free-fall energy (Calvet & Gullbring 1998; Johnston et al. 2011). The very low X-ray luminosity of $<0.01L_{\odot }$ (Anderson et al. 2011) suggests that the free-fall energy is largely absorbed in the stellar surface. Hence, the accretion luminosity that goes into heating the stellar surface is given by

The model stellar luminosity, $L_{bol}$, including the luminosity of the protostar, $L_{\ast }$, and the accretion heating, $L_{heat}$, is

A stellar luminosity of $L_{bol}=1.3\times {10}^4{L}_{\odot }$ determined from infrared observations (Johnston et al. 2011) is applied to constrain the emerging flux at model stellar surface

where $F_{\ast }\equiv L_{\ast }/4\pi R_{\ast }^2$ and is determined by the protostellar mass, $M_{\ast }$, following relations of zero-age main-sequence stars.

The total velocity dispersion of line broadening, $\sigma$, is computed by combining the turbulent broadening of velocity dispersion, ${\sigma }_{nt}$, in quadrature with the thermal broadening, $c_s\equiv \sqrt{kT/2.3m_H}$, where $m_H$ is the mass of Hydrogen atom. For the envelope, we adopted Larson’s law to describe the turbulent broadening of velocity dispersion, ${\sigma }_{nt}=1.10(2r/1pc{)}^{0.38}km{s}^{-1}$ (Larson 1981). Since this turbulent broadening is smaller than our spectral resolution of $0.8km{s}^{-1}$ in most part of the envelope, we apply a line profile with instrumental broadening, ${\Phi }_{\nu }$ (see Appendix A). In the disk, the Shakura-Sunyaev $\alpha$ parameter (Shakura & Sunyaev 1973), is used to describe the kinematic viscosity, ${\nu }_k=\alpha c_{s}H$. In the thin-disk theory, the kinematic viscosity connects mass accretion rate to surface density through

Hence, we have the $\alpha$ parameter described by

which just weakly depends on $T$ and $R$ with $\alpha \propto T^{-1/2}{R}^{-1/4}$. For non-thermal broadening in the disk, we found the ratio of the turbulent to thermal broadening through the $\alpha$ parameter with $\alpha ={\sigma }_{nt}/c_s+B^2/4\pi {\rho }_d{c}_s^2\approx {\sigma }_{nt}/c_s$, where $B$ is the magnetic field. Here we assume the kinematic viscosity mainly attributed to turbulence. The magnetic term becomes important if the field strength is greater than a critical value of $B_{crit}=\sqrt{4\pi {\rho }_d{c}_s^2}\ge 5mG$ at $R_c$. Given the typical field strength of $\le 1mG$ in star-forming cores (Crutcher et al. 2010), the approximation of the $\alpha$ parameter for deriving the turbulent broadening is reasonable in I20126. We kept $\alpha \le 1$ in the disk to make turbulence subsonic when optimizing models.

The model is specified by nine adjustable parameters, listed in Table 2, that are optimized by ${\chi }^2$ minimization. Since our observations are insensitive to $R_{in}$, we do not intend to determine its value but assume $R_{in}=5R_{\ast }$ based on previous infrared studies (Johnston et al. 2011) and findings in disks around low-mass stars (Shu et al. 1994). The best-fit model gives a disk mass of $1.5M_{\odot }$ and a centrifugal radius, $R_c$, of $858AU$ rotating about a $12.0M_{\odot }$ protostar with an accretion luminosity of $2.7\times {10}^3{L}_{\odot }$, giving a disk mass accretion rate, ${\dot{M}}_d$, of $3.9\times {10}^{-5}{M}_{\odot }y{r}^{-1}$. The optimization obtains a reduced ${\overline{\chi }}^2$ value of 1.6, which is optimized with $661696$ pixels including both the continuum and spectral line data, equivalent to $7300$ independent data points. The derived parameters along with their uncertainties are listed in Table 3. The synthetic continuum image from the best-fit model is shown in Fig. 1b, while the integrated intensity images and the P-V diagrams of the $C{H}_{3}CN$ and $C{H}_{3}OH$ emissions are shown in Fig. 3f–3j and Fig. 4f–4j, respectively. For choices of $R_{in}$ in the range of $2R_{\ast }-10R_{\ast }$, parameters of the best-fit models vary slightly within 4% of those listed in Table 2 with the exception of $B_T$, which compensates the change of $R_{in}$ and varies within 30%.

The accretion timescale is estimated by

which gives the time for a mass element in the disk to be accreted, and ${\tau }_{acc}\propto R$. The rotation period of a Keplerian disk is

Hence the ratio ${\tau }_{acc}/P_{rot}\propto R^{-1/2}$ and reaches a minimum at $R_c$. As long as ${\tau }_{acc}>P_{rot}$ at $R_c$, the accretion timescale is longer than the rotation period at all radii, and the accretion is mediated by a rotationally supported disk. In I20126, the surface density at $R_c$ is $2.8gc{m}^{-2}$, which implies an inward velocity, $v_R$, of $0.11km{s}^{-1}$. The accretion timescale at $R_c$, ${\tau }_{acc}=R_c/v_R$, is $3.7\times {10}^{4}yr$, longer than the rotation period of $6.9\times {10}^{3}yr$, so the accretion is mediated by a rotationally supported disk.

As the stellar luminosity is actually a proxy for the stellar mass, our model derives a larger stellar mass, leading to a moderately inclined disk at an angle of $i={48}^{\circ }$ rather than an edge-on geometry, i.e. $i={90}^{\circ }$ (Cesaroni et al. 2005, 2014). A smaller inclination angle of ${41}^{\circ }$ has also been suggested by numerical simulations with a misaligned magnetic field with respect to the disk rotation axis (Shinnaga et al. 2012). The temperature of the disk is above $90K$ at all radii and becomes warmer than the enveloping accretion flow beyond $19AU$. Figure 5 shows the temperature distributions of the disk and the envelope. The mid-plane density of the disk is everywhere larger than $1.6\times {10}^{8}c{m}^{-3}$ and is significantly higher than the critical density for collisional de-excitation for thermalization of the $C{H}_{3}CN$ lines, which is about $5\times {10}^{6}c{m}^{-3}$ at $100K$. The densities in the envelope are lower and thus may not fully satisfy the LTE conditions. Since the observed lines are dominated by the disk component, the approximate treatment of the envelope is not significant.

The derived gas density and mass depend on the molecular abundances, which are assumed to be constant in the entire system, including the disk and envelope. While variations of the molecular abundance cannot be ruled out, the fact that the derived temperature exceeds ice mantle sublimation point readily implies the enhancement of $C{H}_{3}CN$ in the gas phase. The derived $C{H}_{3}CN$ abundance of $2.4\times {10}^{-8}$ is comparable with those found in hot molecular cores (Hernández-Hernández et al. 2014). An abrupt jump of $C{H}_{3}CN$ abundance within the domain of interests is therefore not expected.

Having determined the properties of the disk by comparison with our observations, we can assess its dynamical stability. We calculate the Toomre-$Q$ parameter, $Q=c_s\Omega /\pi G\Sigma$, and find it larger than $2.8$ everywhere in the disk, which makes the disk stable to fragmentation (Fig. 5). The fractional uncertainty of Toomre-$Q$ is about 27% through the disk. The Shakura-Sunyaev $\alpha$ parameter (Shakura & Sunyaev 1973) is $\sim 1$ within $10AU$ and decreases slowly to $0.88$ at $R_c$. The value of $\alpha$ is much higher than the typical value of $\alpha =0.1$, above which fragmentation occurs in simulations of isolated disks with simple cooling laws (Rice et al. 2005). Yet numerical models appropriate for high-mass protostars tend to find higher $\alpha$ (Kratter et al. 2010). A direct comparison between observations and simulations may still be premature because the evolution of the accretion flow in I20126 may not necessarily follow the evolution assumed for the simulations (Kratter et al. 2008).

We present new $345GHz$ continuum and spectral line observations of the disk around the high-mass protostar I20126 with the SMA that achieved the highest angular resolution ($0\text{farcs}37$, equivalent to $\sim 600AU$) to resolve the accretion flow. The continuum emission shows a fairly symmetric morphology, which suggests the disk at a moderate inclination angle. Observations of $C{H}_{3}CN$ and $C{H}_{3}OH$ lines resolve the central $2000AU$ region, where the kinematics display a clear rotation velocity pattern. Position-velocity diagrams of the $C{H}_{3}CN$ lines reveal progressively steeper velocity gradients with increasing upper state energy, indicating an increase in the rotational velocity of the hotter gas nearer the protostar. Such spin-up motions are characteristics of a rotationally supported disk.

We assess the dynamical stability of this massive disk through the Toomre-$Q$ parameter. To evaluate the Toomre-$Q$ as a function of radius through the disk, we measure the gas density, temperature, and rotational velocity in the disk by comparing data with synthetic data generated by radiative transfer models analogous to the standard model of star formation with disk-mediated accretion. Given a luminosity of $1.3\times {10}^4{L}_{\odot }$, the optimized model finds a disk mass of $1.5M_{\odot }$ and a centrifugal radius of $858AU$ rotating about a $12.0M_{\odot }$ protostar with a disk mass accretion rate of $3.9\times {10}^{-5}{M}_{\odot }y{r}^{-1}$. These physical conditions render $Q>2.8$ everywhere in the disk, which makes the disk stable to fragmentation.

Our high angular resolution SMA observations of I20126 provide evidence for a stable massive accretion disk around a high-mass protostar. In contrast to some theoretical expectations of massive disks prone to local instabilities, the disk of I20126 is found to be hot and stable to fragmentation even as the accretion proceeds at a high rate. Such conditions may help to maintain the disk around massive stars and preserve opportunities for developing companions or a planetary system in a later phase of the protostellar evolution.