Grand-design Spiral Arms in a Young Forming Circumstellar Disk

We perform a long-term MHD simulation using a three-dimensional nested-grid code. We refer readers to Machida & Hosokawa (2013) and Machida et al. (2014, 2016) for details. This code solves the MHD equations with self-gravity and Ohmic dissipation. Instead of solving expensive radiation transfer, we adopt the barotropic approximation as in Machida & Matsumoto (2012) to mimic the realistic thermodynamics including the initial isothermal collapse and the quasi-adiabatic first core phase (e.g., Masunaga & Inutsuka 2000; Tomida et al. 2010b, 2013):

$$\begin{eqnarray}&&T=10\,{\rm{K}}\left[1+{\left(\displaystyle \frac{\rho }{{\rho }_{\mathrm{crit}}}\right)}^{2/3}\right],\end{eqnarray} \tag{ 1 }$$

where T and ρ are the gas temperature and density, and ρ_crit = 3.84 × 10⁻¹⁴ g cm⁻³ is the critical density where the gas becomes adiabatic. The ionization degree X_e and Ohmic resistivity η are calculated using simple formulae (Nakano et al. 2002; Machida et al. 2014):

$$\begin{eqnarray}&&\eta =\displaystyle \frac{740}{{X}_{e}}\sqrt{\displaystyle \frac{T}{10\,{\rm{K}}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{X}_{e}=5.7\times {10}^{-4}{\left(\displaystyle \frac{n}{{\mathrm{cm}}^{-3}}\right)}^{-1},\end{eqnarray} \tag{ 3 }$$

where n denotes the gas number density. These are fitting formulae assuming the Mathis–Rumpl–Nordsieck dust size distribution (Mathis et al. 1977) and the ionization rate is dominated by cosmic rays, ζ = 10⁻¹⁷ s⁻¹. This resistivity is conservative (i.e., Ohmic dissipation is weak) as shielding of cosmic rays is not included.

The initial condition is a supercritical Bonnor–Ebert-like sphere (Ebert 1955; Bonnor 1956), whose central density, temperature, radius, and mass are ρ_c = 2.2 × 10⁻¹⁸ g cm⁻³, T = 10 K, R = 6.1 × 10³ au, and M₀ = 1.25 M_⊙, respectively. The sphere has a solid-body rotation around the z-axis with an angular velocity of Ω = 1.5 × 10⁻¹³ s⁻¹ and uniform magnetic field aligned to the rotation axis with B_z = 36 μG. The mass-to-flux ratio averaged over the sphere is μ/μ_crit ∼ 3. These parameters are motivated by observed star-forming clouds (Goodman et al. 1993; Crutcher 2012; Li et al. 2014), but not fine-tuned to reproduce Elias 2–27. The results do not change qualitatively even when we use different parameters.

Each nested-grid level consists of $({N}_{x},{N}_{y},{N}_{z})=(64,64,32)$ cells, with mirror symmetry imposed on the z = 0 plane, and a finer level is created to resolve the local Jeans length with at least eight cells (Truelove et al. 1997). The finest resolution achieved around the center of the computational domain is 0.75 au, while typical resolution at the disk scale (R ∼ 100–200 au) is about 3–6 au. In order to avoid prohibitively small time steps after formation of a protostellar core, we insert a sink particle with an accretion radius of 1 au at the center of the domain (Machida et al. 2014).

We run the simulation as long as possible, and at the end of the simulation the protostar age reaches about 47,000 years. The masses of different components are shown in Figure 1(a). Here, we define the outflow (M_outflow) to be where the gas has a radial velocity larger than 10% of the local isothermal sound speed (${v}_{r}\gt 0.1{c}_{{\rm{s}},\mathrm{iso}}$). The disk (M_disk) is identified as the gas exceeding the critical density ρ_d, which is defined as the minimum density within a region where the rotational motion dominates the radial motion (${v}_{\phi }\gt 2| {v}_{r}|$ and v_ϕ > 0.6 v_Kep where v_Kep is the local Keplerian speed) and not vertically outflowing (${v}_{z}\lt 0.1{c}_{{\rm{s}},\mathrm{iso}}$). The protostar mass M_* is the mass absorbed by the sink particle. The (bound) envelope mass is the total of the rest calculated within the initial cloud radius; ${M}_{\mathrm{env}}={M}_{0}-{M}_{* }-{M}_{\mathrm{disk}}-{M}_{\mathrm{outflow}}$. The object is almost near the end of the Class-I phase as the envelope mass gets as low as 15% of the initial cloud mass. The disk-to-star mass ratio remains almost constant and is about 30%–40% throughout the Class-0 and I phases. Assuming that all the disk material will eventually accrete onto the star, the star formation efficiency will be ∼50%.

Zoom In Zoom Out Reset image size

In the earliest phase, the disk size remains small because magnetic fields transport angular momentum efficiently. However, as the disk evolves, the magnetic angular momentum transport becomes less dominant because magnetic fields dissipate and also because the envelope density decreases and magnetic braking becomes less efficient. Then the disk becomes gravitationally unstable, and grand-design m = 2 spiral arms form spontaneously (see also Hennebelle et al. 2016). The gravitational torque between these non-axisymmetric structures transports angular momentum efficiently and controls the disk evolution. As often pointed out, these material spiral arms wind up tightly and disappear in several orbits by the shearing rotation. However, the disk becomes gravitationally unstable again and spiral arms form recurrently (see the animation of Figure 2). While the spiral arms exist and are transferring the angular momentum, the disk radius increases and the disk becomes highly eccentric. As they transport angular momentum, the disk stabilizes and circularizes. As a result, the disk radius oscillates while it grows, as seen in Figure 1(b). Toward the end of the Class-I phase, the disk size becomes larger than 200 au in the expanding phases with spiral arms.

In order to quantify this behavior, we plot the Q parameter (Toomre 1964) averaged over the disk, which indicates the disk instability, and the amplitude of the spiral arms in Figure 1(c). Here, the Q parameter is defined as the mass-weighted average over the disk:

$$\begin{eqnarray}&&\langle Q\rangle =\displaystyle \frac{{\int }_{\rho \gt {\rho }_{{\rm{d}}}}\tfrac{{c}_{s}\kappa }{\pi G{\rm{\Sigma }}}{\rm{\Sigma }}{dS}}{{\int }_{\rho \gt {\rho }_{{\rm{d}}}}{\rm{\Sigma }}{dS}},\end{eqnarray} \tag{ 4 }$$

where c_s denotes the local sound speed, κ the epicyclic frequency, G the gravitational constant, and Σ the disk surface density, respectively. We use κ = Ω_Kep assuming the disk rotation is almost Keplerian. The amplitude of the spiral arms is measured by the ratio between the maximum and minimum surface densities integrated over −50 au < z < 50 au within the annulus of 0.70 R_disk < R < 0.75 R_disk excluding the outflowing gas with ${v}_{z}\gt 0.1{c}_{{\rm{s}},\mathrm{iso}}$. The Q value increases gradually, indicating that the disk is stabilized overall as it evolves. However, the outer disk still becomes unstable and spiral arms form repeatedly. The disk stabilizes again as the spiral arms transfer the angular momentum, which is observed as a strong correlation between the spiral arm amplitude and the increase of the Q value.

Figure 2 shows the typical evolution of the disk, emergence and decay of the spiral arms. The disk swings between the two states in about 1500 years, the expanding state with spiral arms and the relatively circular, flat state. The orbital timescale is ${t}_{\mathrm{rot}}\sim {\left(\tfrac{4{\pi }^{2}{R}_{\mathrm{disk}}^{3}}{{{GM}}_{* }}\right)}^{1/2}\sim 1500$ years when M_* = 0.44 M_⊙ and R_disk = 100 au, so the transition timescale corresponds to the disk dynamical timescale. We also show the distribution of the plasma beta and the local Q value in Figure 2. The plasma beta is significantly higher than unity within the disk, indicating the disk is weakly magnetized and angular momentum transport by magnetic fields is not efficient anymore. Spiral arms exhibit low Q (Q ≲ 1), but as the disk expands, the Q value gets higher as the disk is stabilized by the angular momentum transport. Figure 3 shows the azimuthally averaged radial profiles of the dust surface density, temperature, and Q value. Our model has a higher surface density and temperature compared to the estimates of Pérez et al. (2016), but it should be noted that the assumption in Pérez et al. (2016) that the surface density should follow a power-law distribution with a cutoff is strong, and the disk does not necessarily have such a distribution when the disk is unstable and dynamically evolving. The Q value is as low as unity where the spiral arms are present, implying that the disk structure is regulated by the balance between the gravitational instability and the angular momentum transport by the spiral arms. As seen in Figure 1, the recurrent spiral arms persist for a significant fraction of the time until the end of the Class-I phase. Assuming that the spiral arms are visible when the amplitude is larger than 50, the occurrence probability of the spiral arms is about 50%. Therefore, there is sufficiently high probability that we can observe such grand-design material spiral arms if the disk is massive. Observations of Elias 2–27 suggest that its disk can be indeed massive, M_disk ∼ 0.14 M_⊙ (Andrews et al. 2009; Pérez et al. 2016).

Zoom In Zoom Out Reset image size

We calculate evolution of the central protostar using the STELLAR stellar evolution code (Yorke & Bodenheimer 2008; Hosokawa et al. 2013; Sakurai et al. 2015), as in Machida & Hosokawa (2013). Time evolution of the intrinsic luminosity and stellar radius as well as the accretion rate onto the sink particle are shown in Figure 1(d). At the epoch we choose, the mass, radius, and intrinsic luminosity of the star are M_* = 0.444 M_⊙, R_* = 2.935 R_⊙, and L_* = 1.604 L_⊙, respectively. The accretion rate in the MHD simulation is still high, a few ×10⁻⁶ M_⊙ yr⁻¹, and the corresponding accretion luminosity should be as high as L_acc ∼ 14 L_⊙ if the gas accretion occurs steadily. However, this accretion rate is time-averaged and there can be short-time accretion variability due to small-scale unresolved structures. Both theoretically and observationally, it is known that most of young stars are fainter than theoretical estimates assuming the steady accretion rate (the so-called luminosity problem), and this is often interpreted as a result of episodic accretion. That is, while the time-averaged accretion rate is still high, the actual accretion onto the protostar happens in long quiescent phases and short burst phases triggered by some physical instabilities within the disk (e.g., Zhu et al. 2010; Dunham & Vorobyov 2012; Stamatellos et al. 2012). Considering this, we adopt $\dot{M}=8\times {10}^{-8}\,{M}_{\odot }{\mathrm{yr}}^{-1}$ to calculate the accretion luminosity based on the observation of Elias 2–27 (Najita et al. 2015), which gives L_acc = 0.379 L_⊙. Assuming a blackbody spectrum, the effective stellar temperature is T_eff = 4000 K. This agrees well with the observed spectral type of Elias 2–27, which is M0 (Luhman & Rieke 1999; Andrews et al. 2009; Najita et al. 2015).

While the barotropic approximation we adopt in the MHD simulation is convenient, it does not fully take into account the radiation heating and cooling after formation of the protostar. To simulate the observation better, we recalculate the dust temperature distribution assuming radiative equilibrium under irradiation from the central star. For this purpose, we utilize RADMC-3D⁶ (Dullemond 2012), which computes radiation transfer using the Monte Carlo method. In addition to the stellar radiation, we impose an external interstellar radiation field of $T=10\,{\rm{K}}$.

Before the radiation transfer simulation, we modify the density distribution. First, we remove the gas and dust outside the initial cloud. Second, because the sink particle only absorbs the gas within the accretion radius whose density is higher than a critical density, ρ_sink = 10⁻¹² g cm⁻³, the gas density around the sink remains unphysically high. Because the gas near the sink should already have accreted, we lower the gas density within a cylinder of R = 1.5 au and z = 4.5 au around the sink to ρ = 10⁻¹⁸ g cm⁻³. This removal of the gas affects the temperature near the protostar, but the removed gas mass is only 5 × 10⁻⁵ M_⊙ and the disk temperature in the large scale (R > 10 au) is not sensitive to how we remove the gas.

For dust opacities, we adopt the monochromatic opacity tables of Semenov et al. (2003),⁷ using the composite aggregate dust model of the normal abundance. The dust-to-gas fraction of this model is 0.014 in the cold (T ≲ 150 K) region where all the dust components exist. This model consists of five tables for different temperatures regarding dust evaporation. To take the dust evaporation into account, we run RADMC-3D repeatedly. First, we run RADMC-3D using the opacity table for the lowest temperature everywhere. Based on the result, we introduce the next opacity table for the higher temperature where the temperature is higher than the evaporation temperature. We repeat this procedure until distributions of all the dust components considering their evaporation are obtained. Then we run RADMC-3D again using 2 × 10⁹ photons to make the temperature distribution as smooth as possible. We include anisotropic scattering and do not use the modified random walk algorithm. Note that our opacity and resistivity are not fully consistent because the dust evaporation is not considered in the resistivity.

The resulting temperature distribution is shown in Figure 4, and the azimuthally averaged radial profile is shown in Figure 3(b). The resulting disk temperature in the outer region (R ≳ 100 au) is about 10–30 K, which is close to the temperature often used in observational studies. The temperature obtained from the radiation transfer calculation is lower than that of the barotropic approximation within the disk near the central star, while it is higher in the envelope. Although the temperature affects the disk stability, we expect that the overall disk evolution, especially in the outer region, is simulated reasonably well because the typical temperature in the outer disk is 10–20 K, which is close to the temperature from the barotropic approximation.

Zoom In Zoom Out Reset image size

The temperatures obtained from the barotropic approximation and radiation transfer calculation are considerably different, but we still can expect the gravitational instability if we properly take radiation transfer into account in our MHD simulation because the temperature from the radiation transfer calculation is lower than the barotropic approximation in most regions within the disk. Also, the Q value has a dependency on the temperature as Q ∝ T^1/2 and therefore the difference in Q should be within a factor of 2. It should be noted, however, we ultimately need radiation MHD simulations to model the disk evolution accurately.

For comparison with the ALMA observation of Elias 2–27, we calculate the intensity distribution of dust continuum at 1.3 mm again using RADMC-3D. The opacity table we use has ${\kappa }_{1.3\mathrm{mm}}\sim 0.95\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ per unit dust mass in most regions within the disk. Note that this value is lower than the opacity used in Pérez et al. (2016) by a factor of ∼2.5. To match the observation, we set the inclination angle to be 558. In order to avoid unphysical noise in the ALMA simulation due to the edges of the image, the image size must be larger than the primary beam size of ALMA, which is about 27 arcsec at 1.3 mm. Assuming the object is located at a distance of 139 pc, we calculate the image covering (10,000 au)² with 4000² pixels. The result is shown in the left panel of Figure 5.

Zoom In Zoom Out Reset image size

The column density along the line of the sight from the edge to the center of the cloud, excluding the gas within the sink radius, is ${N}_{{{\rm{H}}}_{2}}\sim 6.34\times {10}^{21}\,{\mathrm{cm}}^{-2}$, indicating that the envelope is almost exhausted. This value is lower than one inferred from the observed visual extinction A_V ∼ 14.8 (Andrews et al. 2009; Najita et al. 2015). However, the column density estimate in our simulation is only a lower limit because we exclude the structure near the sink and outside the initial cloud.

Then we perform simulated observation using the simobserve and simanalyze tasks of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007). We chose observation parameters as close to Pérez et al. (2016) as possible; the bandwidth is 6.8 GHz, the integration time is 12.5 minutes, and the source position is R.A. (J2000) = 16^h26^m45024, decl. (J2000) = −24^d23^m08250. We use the "alma.cycle4.5.cfg" array configuration file and the Briggs waiting with the robustness of 0.5. This produces a synthetic beam of 0.29 arcsec × 0.26 arcsec, which is very close to the original observation. The result of the synthetic observation is shown in Figure 5, as well as the original observation of Elias 2–27 (Pérez et al. 2016).

The disk structure agrees well with the ALMA observation of Elias 2–27. The central bright disk and the grand-design spiral arms are clearly visible. The disk size, brightness, and thickness of the spiral arms are largely consistent with the observation. Therefore, we conclude that the observed spiral arms can be well explained by material spiral arms formed by gravitational instability.