Dynamical Gaseous Rings in Global Simulations of Protoplanetary Disk Formation

The equations of mass, momentum, and energy transport in the thin-disk limit are

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}=-{{\boldsymbol{\nabla }}}_{p}\cdot \left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\right)+{\left[{\boldsymbol{\nabla }}\cdot \left({\rm{\Sigma }}{{\boldsymbol{v}}}_{p}\otimes {{\boldsymbol{v}}}_{p}\right)\right]}_{p}\\ \quad \ =\ -{{\boldsymbol{\nabla }}}_{p}P+{\rm{\Sigma }}{{\boldsymbol{g}}}_{p}+{\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{\Pi }}\right)}_{p},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial e}{\partial t}+{{\boldsymbol{\nabla }}}_{{\rm{p}}}\cdot \left({{ev}}_{p}\right)\\ \quad \ =\ -P\left({{\boldsymbol{\nabla }}}_{p}\cdot {v}_{p}\right)-{\rm{\Lambda }}+{\rm{\Gamma }}+{\left({\boldsymbol{\nabla }}\cdot v\right)}_{{pp}^{\prime} }:{{\boldsymbol{\Pi }}}_{{pp}^{\prime} },\end{array}\end{eqnarray} \tag{ 3 }$$

where subscripts p and p' refer to the planar components (r, ϕ) in polar coordinates, Σ is the mass surface density, e is the internal energy per surface area, P is the vertically integrated gas pressure calculated via the ideal equation of state as P = (γ − 1)e with γ = 7/5, ${{\boldsymbol{v}}}_{p}={v}_{{\rm{r}}}\hat{{\boldsymbol{r}}}+{v}_{\phi }\hat{{\boldsymbol{\phi }}}$ is the velocity in the disk plane, ${{\boldsymbol{g}}}_{p}={g}_{{\rm{r}}}\hat{{\boldsymbol{r}}}+{g}_{\phi }\hat{{\boldsymbol{\phi }}}$ is the gravitational acceleration in the disk plane, and ${{\boldsymbol{\nabla }}}_{{\rm{p}}}=\hat{{\boldsymbol{r}}}\partial /\partial r+\hat{{\boldsymbol{\phi }}}{r}^{-1}\partial /\partial \phi$ is the gradient along the planar coordinates of the disk.

Turbulent viscosity is taken into account via the viscous stress tensor

$$\begin{eqnarray}&&{\boldsymbol{\Pi }}=2{\rm{\Sigma }}\nu \left({\rm{\nabla }}v-\displaystyle \frac{1}{3}({\rm{\nabla }}\cdot v){\boldsymbol{e}}\right),\end{eqnarray} \tag{ 4 }$$

where ν is the kinematic viscosity, ∇v is a symmetrized velocity gradient tensor, and ${\boldsymbol{e}}$ is a unit tensor. We parameterize the magnitude of kinematic viscosity ν = α_effc_sH, where c_s is the sound speed and H is the vertical scale height, using the α-prescription of Shakura & Sunyaev (1973) with an effective and adaptive α-parameter expressed in the following form to take the layered disk structure into account (Bae et al. 2014):

$$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{a}}}{\alpha }_{{\rm{a}}}+{{\rm{\Sigma }}}_{{\rm{d}}}{\alpha }_{{\rm{d}}}}{{\rm{\Sigma }}},\end{eqnarray} \tag{ 5 }$$

where Σ_a is the gas surface density of the upper MRI-active layer of the disk and Σ_d is the gas surface density of the lower MRI-dead layer. The total surface density of gas is then Σ = Σ_a + Σ_d. We note that Σ_d is set equal to zero if Σ ≤ Σ_a. The parameters α_a and α_d represent the strength of turbulent viscosity in the MRI-active and MRI-dead layers of the disk. Following Bae et al. (2014), we set α_a = 0.01. Although the actual values may vary in wide limits, the averaged value is unlikely to significantly exceed this limit (e.g., Yang et al. 2018). Indeed, numerical hydrodynamics simulations of the disk evolution with a constant α_a = 0.1 yielded disk lifetimes much smaller than the observationally inferred values of 2–3 Myr (Vorobyov & Basu 2009). Furthermore, we define α_d as

$$\begin{eqnarray}&&{\alpha }_{{\rm{d}}}={\alpha }_{\mathrm{MRI},{\rm{d}}}+{\alpha }_{\mathrm{rd}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}{\alpha }_{\mathrm{MRI},{\rm{d}}}=\left\{\begin{array}{cc}{\alpha }_{{\rm{a}}}, & \mathrm{if}\,{T}_{\mathrm{mp}}\gt {T}_{\mathrm{crit}}\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 7 }$$

T_crit is the MRI activation temperature due to thermal effects, such as dust sublimation and/or ionization of alkaline metals, and ${T}_{\mathrm{mp}}=\mu { \mathcal P }/({ \mathcal R }$ Σ) is the disk midplane temperature. Here μ = 2.33 is the mean molecular weight and ${ \mathcal R }$ is the universal gas constant. We note that this definition of T_mp neglects possible vertical variations in the gas temperature. For a more accurate calculation of the midplane temperature one needs to solve the radiative transfer equations, as was done in the recent modification of the thin-disk limit by Vorobyov & Pavlyuchenkov (2017).

As was noted in Bae et al. (2014) based on magnetohydrodynamics simulations of Okuzumi & Hirose (2011), the disk dead layer can have some nonzero residual viscosity, with α_rd as large as 10⁻³–10⁻⁵, due to hydrodynamic turbulence driven by the Maxwell stress in the disk active layer. Therefore, we define the residual viscosity as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{rd}}=\min \left({10}^{-5},{\alpha }_{{\rm{a}}}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{a}}}}{{{\rm{\Sigma }}}_{{\rm{d}}}}\right).\end{eqnarray} \tag{ 8 }$$

This expression takes into account the fact that accretion in the dead zone cannot exceed that of the active zone if the nonzero α_rd is due to turbulence propagating from the active layer down to the disk midplane.

The cooling and heating rates Λ and Γ take the disk cooling and heating due to stellar and background irradiation into account based on the analytical solution of the radiation transfer equations in the vertical direction (see Dong et al. 2016, for details):⁵

$$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \frac{8{\tau }_{{\rm{P}}}\sigma {T}_{\mathrm{mp}}^{4}}{1+2{\tau }_{{\rm{P}}}+\tfrac{3}{2}{\tau }_{{\rm{R}}}{\tau }_{{\rm{P}}}},\end{eqnarray} \tag{ 9 }$$

where σ is the Stefan–Boltzmann constant, τ_R and τ_P are the Rosseland and Planck optical depths to the disk midplane, and κ_P and κ_R (in cm² g⁻¹) are the Planck and Rosseland mean opacities taken from Semenov et al. (2003).

The heating function per surface area of the disk is expressed as

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{8{\tau }_{{\rm{P}}}\sigma {T}_{\mathrm{irr}}^{4}}{1+2{\tau }_{{\rm{P}}}+\tfrac{3}{2}{\tau }_{{\rm{R}}}{\tau }_{{\rm{P}}}},\end{eqnarray} \tag{ 10 }$$

where T_irr is the irradiation temperature at the disk surface determined from the stellar and background blackbody irradiation as

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}^{4}={T}_{\mathrm{bg}}^{4}+\displaystyle \frac{{F}_{\mathrm{irr}}(r)}{\sigma },\end{eqnarray} \tag{ 11 }$$

where ${F}_{\mathrm{irr}}(r)$ is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central star. The latter quantity is calculated as

$$\begin{eqnarray}&&{F}_{\mathrm{irr}}(r)=\displaystyle \frac{{L}_{* }}{4\pi {r}^{2}}\cos {\gamma }_{\mathrm{irr}},\end{eqnarray} \tag{ 12 }$$

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r. The incidence angle is calculated using a flaring disk surface as described in Vorobyov & Basu (2010). The stellar luminosity L_* is the sum of the accretion luminosity ${L}_{* ,\mathrm{accr}}=(1-\epsilon ){{GM}}_{* }\dot{M}/2{R}_{* }$ arising from the gravitational energy of accreted gas and the photospheric luminosity L_*,ph due to gravitational compression and deuterium burning in the stellar interior. The stellar mass M_* and accretion rate onto the star $\dot{M}$ are determined using the amount of gas passing through the inner computational boundary. The properties of the forming protostar (L_*,ph and radius R_*) are calculated using the stellar evolution tracks provided by the software STELLAR (Yorke & Bodenheimer 2008) code. The fraction of accretion energy absorbed by the star is set to 0.05.

The choice of the inner boundary condition in disk formation simulations always presents a certain challenge. Ideally, the inner boundary should be placed near the stellar surface, where the inner disk joins the stellar magnetosphere. Such a small inner boundary is, however, computationally prohibitive even for 2D time-explicit hydrodynamics simulations, as our own, because of strict time step limitations imposed by the Courant condition (Courant et al. 1928) in Keplerian disks. This necessitates the use of a sink cell that cuts out the inner disk region, thus moving the inner disk boundary to a larger distance and relaxing the Courant condition. The farther the sink–disk interface from the star, the longer time step the code can manage and the longer integration times (or higher numerical resolution) one can afford. This, however, comes at a price—too large a sink cell may cut out the disk regions that are dynamically important for the entire disk evolution. In this work, the inner boundary of the computational domain was placed at r_sc = 0.4 au, which allowed us to capture the disk behavior occurring at sub-astronomical-unit scale. To the best of our knowledge, this is the smallest sink cell used in global collapse simulations of prestellar cores. Moving the inner boundary to several astronomical units, as in most collapse simulations, would not allow us to capture the MRI triggering, which plays a key role in the evolution of the inner disk.

Another complication involves the type of the inner boundary. The inner boundary cannot be of reflecting type because it must allow for matter to flow to the sink cell and ultimately on the nascent protostar. If the inner boundary allows for matter to flow only in one direction, i.e., from the disk to the sink cell, then any wave-like motions near the inner boundary, such as those triggered by spiral density waves in the disk, would result in a disproportionate flow through the sink–disk interface. As a result, an artificial depression in the gas density near the inner boundary develops in the course of time because of the lack of compensating backflow from the sink to the disk.

A solution to this problem is to allow matter to flow freely from the disk to the central sink cell and vice versa. This inflow–outflow boundary condition was presented in Vorobyov et al. (2018) and is schematically illustrated in Figure 1. The mass of material that passes to the sink cell through the sink–disk interface is further redistributed between the central protostar and the sink cell as ${\rm{\Delta }}{M}_{\mathrm{flow}}={\rm{\Delta }}{M}_{* }+{\rm{\Delta }}{M}_{{\rm{s}}.{\rm{c}}.}$ according to the following algorithm (note that ΔM_flow is always positive definite):

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{if}\,{{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}\lt {\overline{{\rm{\Sigma }}}}_{\mathrm{in}.\mathrm{disk}}^{n}\,\mathrm{and}\,{v}_{r}({r}_{{\rm{s}}.{\rm{c}}.})\lt 0\,\mathrm{then}\\ \ \ \ {{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}}^{n+1}={{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}+{\rm{\Delta }}{M}_{{\rm{s}}.{\rm{c}}.}/{S}_{{\rm{s}}.{\rm{c}}.}\\ \ \ \ {M}_{* }^{n+1}={M}_{* }^{n}+{\rm{\Delta }}{M}_{* }\\ \mathrm{if}\,{{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}\lt {\overline{{\rm{\Sigma }}}}_{\mathrm{in}.\mathrm{disk}}^{n}\,\mathrm{and}\,{v}_{r}({r}_{{\rm{s}}.{\rm{c}}.})\geqslant 0\,\mathrm{then}\\ \ \ \ {{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}}^{n+1}={{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}-{\rm{\Delta }}{M}_{\mathrm{flow}}/{S}_{{\rm{s}}.{\rm{c}}.}\\ \ \ \ {M}_{* }^{n+1}={M}_{* }^{n}\\ \mathrm{if}\,{{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}\geqslant {\overline{{\rm{\Sigma }}}}_{\mathrm{in}.\mathrm{disk}}^{n}\,\mathrm{and}\,{v}_{r}({r}_{{\rm{s}}.{\rm{c}}.})\lt 0\,\mathrm{then}\\ \ \ \ {{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n+1}={{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}\\ \ \ \ {M}_{* }^{n+1}={M}_{* }^{n}+{\rm{\Delta }}{M}_{\mathrm{flow}}\\ \mathrm{if}\,{{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}\geqslant {\overline{{\rm{\Sigma }}}}_{\mathrm{in}.\mathrm{disk}}^{n}\,\mathrm{and}\,{v}_{r}({r}_{{\rm{s}}.{\rm{c}}.})\geqslant 0\,\mathrm{then}\\ \ \ \ {{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n+1}={{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n}-{\rm{\Delta }}{M}_{\mathrm{flow}}/{S}_{{\rm{s}}.{\rm{c}}.}\\ \ \ \ {M}_{* }^{n+1}={M}_{* }^{n}.\end{array}\end{eqnarray*}$$

Here Σ_s.c. is the surface density of gas in the sink cell, ${\overline{{\rm{\Sigma }}}}_{\mathrm{in}.\mathrm{disk}}$ is the averaged surface density of gas in the inner active disk (the averaging is usually done over 1 au immediately adjacent to the sink cell), S_s.c. is the surface area of the sink cell, and v_r(r_s.c.) is the radial component of velocity at the sink–disk interface. We note that v_r(r_s.c.) < 0 when the gas flows from the active disk to the sink cell and v_r(r_s.c.) > 0 in the opposite case. The superscripts n and n + 1 denote the current and the updated (next time step) quantities. The exact partition between ΔM_* and ΔM_s.c. is usually set to 95%:5%. The influence of the ΔM_*:ΔM_s.c. partition on the disk evolution is studied in Vorobyov et al. (2019).

Zoom In Zoom Out Reset image size

The calculated surface densities in the sink cell ${{\rm{\Sigma }}}_{{\rm{s}}.{\rm{c}}.}^{n+1}$ are used at the next time step as the inner boundary values for the surface density. The radial velocity and internal energy at the inner boundary are determined from the zero gradient condition, while the azimuthal velocity is extrapolated from the active disk to the sink cell assuming a Keplerian rotation. These inflow–outflow boundary conditions enable a smooth transition of the surface density and angular momentum between the inner active disk and the sink cell, preventing (or greatly reducing) the formation of an artificial drop in the surface density near the inner boundary. We ensure that our boundary conditions conserve the total mass budget in the system. Finally, we note that the outer boundary condition is set to a standard free outflow, allowing for material to flow out of the computational domain, but not allowing any material to flow in.

Our numerical simulations start from the gravitational collapse of a starless cloud core and continue into the embedded phase of star formation when a star, disk, and envelope are formed. The initial surface density and angular velocity of the cloud core are similar to those derived by Basu (1997), for axially symmetric core collapse

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{{r}_{0}{{\rm{\Sigma }}}_{0}}{\sqrt{{r}^{2}+{r}_{0}^{2}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rm{\Omega }}=2{{\rm{\Omega }}}_{0}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\left[\sqrt{1+{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{2}}-1\right].\end{eqnarray} \tag{ 14 }$$

Here Σ₀ and Ω₀ are the corresponding values at the center of the core and r₀ is the radius of the central plateau. These parameters depend on the intended initial mass and rotational energy of the core (Vorobyov & Basu 2010). The simulations are terminated after about 0.5 Myr in the T Tauri phase, when most of the envelope has accreted onto the central star plus disk system. We conducted a total of 10 disk simulations, which are listed in Table 1. The prescript model1 or model2 corresponds to the total mass of the gas in the initial prestellar cloud core, which is set at 1.152 and 0.346 M_⊙, respectively. With the two parent core masses we aim to demonstrate the dependence of the disk structure and evolution on the final stellar mass. With a larger core mass model1 should ultimately form a Sun-like star and model2, which starts with much smaller cloud core mass, should form a fully convective, sub-solar-mass star. All of the initial cloud cores were constructed such that the ratio of the rotational to the gravitational energy, β, was approximately 0.135%. This value is consistent with the observations of prestellar cores (Caselli et al. 2002).

The MRI-driven turbulence is thought to be one of the dominant sources of angular momentum transport in protoplanetary disks (Hawley et al. 1995; Turner et al. 2014). This mechanism acts in a shearing Keplerian disk when the gas in the disk is sufficiently ionized and weakly coupled to the magnetic field. The midplane temperature in a typical protoplanetary disk is often not high enough to sustain collisional ionization, and thus it forms a magnetically dead zone. The critical temperature, T_crit, corresponds to the threshold temperature at which the gas in the disk is sufficiently ionized, so that the disk makes an effective transition from an MRI-dead state to an MRI-active one. The exact value of T_crit is difficult to calculate, as it depends on local conditions in the disk such as the distribution of alkali metals, ionization rate, abundance of small grains, and threshold ionization fraction for MRI turbulence. However, there is a general consensus that for protoplanetary disks the transition occurs abruptly at ⪆1000 K when there is an almost exponential rise in the electron fraction in the gas, primarily due to the ionization of potassium (Umebayashi 1983; Gammie 1996). In accordance with recent investigations (e.g., Zhu et al. 2010; Bae et al. 2014), we conducted simulations with two values of T_crit, 1300 and 1500 K, which is denoted in the model name by the number after letter "T."

Another uncertainty in the model arises owing to the thickness of the active layer, which depends on the level of disk ionization. The cosmic rays are considered to be the dominant source of ionization, penetrating an approximately constant column density of about 100 g cm⁻² throughout the disk (Umebayashi & Nakano 1981), which is a canonical value for layered disk models. Hard X-rays (∼5 keV) originating from the star are usually neglected, since their penetration depth is typically an order of magnitude smaller (Igea & Glassgold 1999). However, the cosmic rays may not always reach the disk surface unattenuated. Similar to the action of solar wind (Spitzer & Tomasko 1968), it is possible that the powerful stellar winds or disk winds of a young star can shield the inner disk from the high-energy cosmic rays. The magnetic field structure of the system can also hinder the propagation of the cosmic rays (Padovani & Galli 2011). These factors could result in a shielding effect, significantly diminishing active layer thickness. Thus, in our simulations we used two values of Σ_a, 100 and 10 g cm⁻², representing an unshielded and shielded disk, respectively. The value of Σ_a is denoted in the model name by the number after letter "S."

The first model, model1_T1300_S100, is treated as a representative, fiducial magnetically layered disk model, and model1_const_alpha is treated as a fiducial fully MRI-active disk model in this paper. The two models are presented in detail, and all other models are evaluated in comparison with these fiducial models. Apart from the three parameters—M_gas, T_crit, and Σ_a—all other model parameters are kept identical across the 10 simulations.

The radial and azimuthal resolution of the computational grid was 512 × 512 for all simulations, with logarithmic spacing in the radial direction and uniform spacing in the azimuthal direction. The highest grid resolution in the radial direction was 8.241 × 10⁻³ au at the innermost cell of the disk (0.4167 au), while the poorest resolution near the outer boundary of the cloud (10210 au for model1) was about 202.0 au. We assured the sufficiency of the above resolution, as well as the numerical convergence of the simulations, by comparing global disk properties of the fiducial models at a lower resolution. The values of α_a and α_d in the layered disk models were 10⁻² and 10⁻⁵, respectively, while the fully MRI-active disk models were evolved with a single value of α = 10⁻². The initial gas temperature of the prestellar cloud core was set to 15 K, which was also the temperature of the background radiation.

In this section we elaborate on the properties of the gaseous rings formed inside the dead zones in layered disk models. The inner disk structure of the fiducial layered disk model (model1_T1300_S100) is compared with that of the fiducial fully MRI-active disk model (model1_const_alpha). Both models were evolved from identical initial conditions, and the only difference between them was the viscosity prescription as described in Section 2.3. The classical approach for observational modeling of disks is to consider a Shakura–Sunyaev prescription of viscosity with a constant value of α-parameter, which in our case represents a fully MRI-active disk. Although such a prescription has been demonstrated to work very well in certain astrophysical objects (e.g., disks in cataclysmic variables, X-ray binaries; see Frank et al. 2002), a magnetically layered protoplanetary disk is better represented with an adaptive α. The motivation behind this comparison between the two fiducial models was to demonstrate the differences and their implications, especially in the inner disk smaller than about 15 au, which is the region most interesting for planet formation.

Figure 2 shows the global evolution of the gas surface density distribution in the disk for the two fiducial models over the inner 500 × 500 au region. The prestellar cloud core was constructed such that it was gravitationally unstable, and the time was measured from the beginning of the simulations when it began to collapse. The disk was evolved up to about 0.5 Myr in both simulations, when it can be considered as an early T Tauri object. Note that the two models look qualitatively similar over large scale. The initial phases, when the disk is massive, are characterized by gravitational instability (GI) and fragmentation. Both disks show irregular spiral structure and formation of clumps in this phase. As the evolution progressed, the disks grew and viscously spread, ultimately showing a smoother profile. The major difference between the two models is that the central region of the layered disk model showed a consistently higher surface density as compared to the fully active disk. The overall size of the disk is smaller in the layered disk model, indicating less viscous spread as compared to the fully MRI-active model. As a consequence of relatively larger transport of angular momentum to the outer disk region, the fully MRI-active disk also shows more clump formation through fragmentation.

Zoom In Zoom Out Reset image size

In Figure 3 we focus on the inner region of the fiducial disks at 0.35 Myr. This time corresponds to the fourth column in Figure 2, when the disk structures are typical for both simulations. At this point, the disks had started to settle down, became increasingly symmetric, and showed less flocculent spiral structure. Figure 3 compares a snapshot of four quantities of the fiducial disks—total surface density, α_eff as given by Equation (5), midplane temperature, and Toomre's Q-parameter (calculated as Q = c_sΩ/πGΣ). In each row, the two images on the right-hand side show the distribution of the variables in the inner 10 × 10 au box. On the left-hand side in each row the azimuthally averaged profiles of the same variables are plotted for comparison in 100 au radius. The shaded region shows the extent between the maximum and minimum values at a given radius. Note that the abscissas in the azimuthal profiles are in logarithmic units.

Zoom In Zoom Out Reset image size

We define the dead zone in the layered disk model as the region where α_eff falls below 80% of the fully MRI-active value of α_a = 0.01. The physical motivation behind this definition was that below this threshold of α_eff, the gas surface density in all of the layered models started to diverge from a fully MRI-active disk, as can also be observed in Figure 3. It was possible to consider the dead zone in terms of the surface density, as the dead zone exists only when Σ ≥ Σ_a. However, we will show that this was always an overestimate, and the above definition of dead zone based on the viscosity prescription (α_eff) reflected the status of the inner disk more accurately. The vertical black line in the azimuthal profiles in Figure 3 marks the outer extent of the dead zone at this time, which was at about 16 au. The azimuthal profiles of the quantities from the two simulations in Figure 3 tend to converge beyond this point, and within the dead zone they differed significantly.

The most compelling difference between the fiducial models was the formation of concentric, axisymmetric, gaseous "rings" within the dead zone in the layered disk model. As seen in azimuthal plots, as well as the spatial distributions, the rings offer a remarkable contrast in the local conditions as compared to the fully magnetically active disk. These rings are characterized by a much higher surface density and a lower effective viscosity, indicated by a low value of α_eff. The two vertical lines at about 1.2 and 3.2 au in the azimuthal plots indicate the two maxima in the surface density and hence the position of the two rings. For the most prominent ring, the surface density was about two orders of magnitude higher and the local α_eff was about two orders of magnitude lower as compared to the corresponding values in the fully active disk. The outer ring, although clearly distinguishable, was about an order of magnitude less pronounced, in terms of both Σ and α_eff, as compared to the inner ring. The midplane temperature at the locations of the rings was lower, as compared to a fully active disk, because of the lower viscous heating in these regions.

Toomre's Q-parameter compares the destabilizing effects of self-gravity in the disk against the stabilizing effects of pressure and rotation. This parameter is expected to increase as we get closer to the star owing to the increased rotational velocity and rising temperature, as we see in the fully MRI-active model. However, due to the accumulation of gas, as well as the lower midplane temperature, the Q-parameter was low (≈2) in the rings, which indicates that this region was marginally gravitationally unstable.

Although the properties of the inner disk region in Section 3.1 are typical for the layered disk model, the rings themselves are highly dynamical and constantly evolving structures. In order to demonstrate their behavior across time, we plot spacetime diagrams in Figure 4. This figure compares the azimuthally averaged profiles of quantities Σ, α_eff, and T, the minimum in Toomre's Q-parameter, as well as the viscous torque (τ_ν). Each row shows the evolution of the variables for the fiducial layered disk model on the left and the fully MRI-active disk model on the right. Note that the ordinate (R) is in logarithmic units. The azimuthally averaged profiles were printed every 1000 yr in our numerical code, which determines the temporal resolution in this plot. The time spacing of the Q-parameter was 2000 yr, since the minimum value was calculated using 2D output files. We found that this resolution in time was sufficient for analyzing the dynamical behavior of the disk. The vertical dashed line at 0.35 Myr denotes the snapshot in time corresponding to Figure 3. In accordance with the earlier discussion, the region outside the dead zone was similar for both fiducial models throughout the evolution, while the inner region appeared markedly different.

Zoom In Zoom Out Reset image size

Consider the first row in Figure 4, depicting the evolution of surface density. The green contour plotted at Σ = 1 g cm⁻² shows the approximate extent of the disk, and the yellow contour is plotted at Σ_a = 100 g cm⁻². With the fully MRI-active disk, the surface density approximately monotonically decreases with radius, while also decreasing across time as the disk evolves. Some variations in the surface density at larger distances and earlier times are the result of clump formation through disk fragmentation, in both simulations. In the layered disk model, enhancements in surface density in the region between 0.6 and 5 au started occurring as soon as the disk was formed at 0.08 Myr, while axisymmetric structures resembling rings started appearing after about 0.2 Myr. The rings can be easily identified as the diagonal high surface density bands. As the time evolved, the rings became denser and also migrated closer to the star. Multiple rings could form at a given radius; however, the innermost ring was typically the most prominent. After the inward migration, we can observe that a ring terminated with a sudden discontinuity. We found that at these points in time MRI was triggered, and as a result the ring became unstable. The material in the ring was accreted onto the star, and the resulting eruption qualitatively resembled an FUor outburst. As mentioned earlier, in this paper we will concentrate on the structure and evolution of the inner disk, and the detailed analysis of the outbursts will be conducted in a subsequent paper.

The second row of Figure 4 compares the effective α-parameter, which has a constant value of 0.01 for the fully MRI-active disk. The black contour shows the extent of the dead zone in the layered disk model, defined as the region where the α_eff falls below 80% of this maximum value. The yellow contours depict a level of 10⁻⁴. Thus, within the vicinity of the rings, α_eff is less than 1% of its MRI-active value. This low-viscosity region coincides with the rings observed in the gas surface density. Note that for the layered disk model, at early times the dead zone was consistent with the Σ_a = 100 g cm⁻² contour. As the disk evolved, the gas surface density became exceedingly flatter, and the dead zone was contained within a much smaller radius. The third row compares the azimuthally averaged midplane temperature. For the fully MRI-active disk the temperature distribution mirrored the surface density, which was approximately monotonically decreasing with the radius. The disk was initially hot and then cooled as it evolved. In contrast, the innermost regions of the layered disk model (<1 au) remained at approximately constant temperature throughout its 0.5 Myr evolution. The temperature near the rings was lower than the surroundings because of the decreased viscous heating. The yellow contour shows T_crit = 1300 K. Notice that the midplane temperature in the inner disk periodically spiked above this value for a short period in the layered disk model. The increases in the inner disk temperature coincide with the disappearance of the rings and were thus associated with the MRI outbursts.

The fourth row of Figure 4 compares the minimum in the Toomre's Q-parameter for the two fiducial simulations, which indicates whether the disk is gravitationally unstable at a given radius. The yellow shaded region shows the region where Q < 1, while the green contours show the Q = 2 levels. In both simulations, the initial gravitational collapse phase of a prestellar core can be seen before the disk formation at about 0.08 Myr as gravitationally unstable at all radii. As expected, disks in both simulations were gravitationally unstable at large distances (⪆10 au) and early times (⪅0.35 Myr), which points to the formation of clumps via disk fragmentation. In the vicinity of the rings, the yellow pigments suggest that the rings showed brief instances of GI fragmentation, although not all such occurrences were captured owing to the limited time resolution. In this region the Q-parameter remained consistently low owing to the large total surface density, as well as a low temperature. We elaborate on the role of GI and associated spirals on the evolution of the rings later in this section.

The last row of Figure 4 shows the viscous torque for the two fiducial models, calculated as ${\tau }_{\nu }{(r,\phi )=r({\rm{\nabla }}\cdot {\boldsymbol{\Pi }})}_{\phi }S(r,\phi )$, where ${\rm{\nabla }}\cdot {\boldsymbol{\Pi }}$ is the viscous force per unit area (Equation (2)), r is the lever arm, and S(r, ϕ) is the surface area occupied by a cell with coordinates (r, ϕ) (see Vorobyov & Basu 2009 for full expression). The outer regions are noisy owing to the action of clumps and spirals, especially in the early times. Neglecting the initial ~0.15 Myr, the general direction of the viscous torques in the inner disk (⪅10 au) is highlighted by arrows. The fully MRI-active disk shows a smooth negative viscous torque in this region, indicating a steady inward mass transfer. The evolution of τ_ν in the fiducial layered disk model is again much different; the viscous torques change sign across the rings and act toward building up the rings from both radial directions. The action of viscous torques can explain how the rings form in the inner disk. The outer edge of the dead zone usually undergoes a smooth and gradual transition. As the surface density decreases in the radial direction, the shift from dead to fully active region occurs over several tens of astronomical units. Thus, the mass does not accumulate near the outer boundary of the dead zone. The inner boundary of the dead zone, on the other hand, undergoes a sharp transition when the temperature exceeds T_crit (1300 K for the fiducial model). The viscous torque has a component proportional to the negative of the gradient of kinematic viscosity, so the material is accelerated outward at the inner boundary. In the fiducial layered disk model, the viscosity gradient was relatively large in the vicinity of the rings, as can be inferred from the azimuthal plots of α_eff in Figure 3. Thus, when the gas reached the inner regions of the dead zone, the viscous torques worked in the direction of reinforcing the pileup of material. The resulting positive feedback between viscous torques and gradients of kinematic viscosity (due to the increased surface density) generated the stable ring structures. The viscous torque is proportional to the distance from the axis of rotation; hence, its magnitude decreased in general as the rings moved inward.

The interplay between the magnetically layered structure in the form of rings and gravitational instabilities showed interesting, time-dependent manifestations over a large extent of the disk. Due to the action of viscous torques, the rings accumulated enough mass to show occasional gravitational fragmentation. These fragments could not survive for long owing to strong Keplerian shear at these radii, as well as their proximity to other overdense regions of the ring. However, this threshold for GI fragmentation (Q = 1) captures only an extreme case of the nonlinear outcome of axisymmetric instability. A disk that is locally stable according to this criterion may still generate nonaxisymmetric, large-scale spiral waves (Hohl 1971; Polyachenko et al. 1997). Thus, the onset of GI may occur at a threshold as high as Q ≈ 2, and associated spiral waves can extend over a large fraction of the disk. When the disk is massive, the evolution of GI may show episodes of intense spiral activity followed by quiescent phases (Lodato & Rice 2005). Similar behavior was observed in our simulations; spiral density waves emanated from periodic asymmetries developed in the rings, and they were particularly strong when GI fragmentation occurred. We demonstrate this in Figure 5, which compares the disk properties when gravitational fragmentation was absent in the rings (snapshot at 0.350 Myr, corresponding to Figure 3) to when the rings showed GI fragmentation with Q < 1 (snapshot at 0.326 Myr, corresponding to Figure 6). The first row shows the fractional amplitudes, $\delta {\rm{\Sigma }}/\langle {\rm{\Sigma }}\rangle$, where the denominator is averaged in the azimuthal direction. The spiral density waves with low azimuthal mode number (usually m = 2) originating at a locally unstable region extended across the entire dead zone, while such waves were over an order of magnitude weaker at other times when the rings were more axially symmetric. Note that the fractional amplitude map for the fiducial fully MRI-active disk showed spirals as well, which can be inferred from the first row of Figure 3. These spirals did not show a time-dependent behavior except for their progressive weakening as the disk evolved (see Vorobyov & Basu 2009 for more details).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The large-scale spiral waves can contribute toward angular momentum transport across the dead zone and thus enable inward migration of the rings. The second row of Figure 5 quantifies the angular momentum transport due to GI via gravitational α-parameter calculated as a ratio of gravitational stress G_xy and vertically averaged pressure P,

$$\begin{eqnarray}&&{\alpha }_{\mathrm{GI}}=\displaystyle \frac{{G}_{{xy}}}{P}=\displaystyle \frac{1}{4\pi {GP}}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial x}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial y},\end{eqnarray} \tag{ 15 }$$

where Φ is gravitational potential. The azimuthally averaged α_GI can be added to α_eff to give an estimate of total viscosity (and thus the angular momentum transport) at a given radius due to both GI and MRI processes. When azimuthal asymmetries were absent in the rings, the inner disk showed consistently low values of α_GI ≈ 0 (the fiducial fully MRI-active disk showed similarly low values). In presence of asymmetries with GI fragmentation, α_GI was much stronger, reaching about 0.2 owing to strong spiral density waves. This value is over an order of magnitude larger than the largest possible value of α_eff = α_a = 0.01 in the layered disk. Thus, these GI spiral waves were capable of redistributing a significant amount of angular momentum and gravitational energy across the dead zone. Note that even when the Q-criterion is fulfilled locally, the disk needs to cool efficiently over a dynamical timescale to form GI fragments (Rafikov 2005). As rings migrated inward, they may be more stable against fragmentation because of inefficient cooling, thus showing a slower rate of migration.

An additional aspect of dynamical behavior of the inner disk was occasional formation of vortices in the vicinity of the rings. Figure 6 depicts an instance of a vortex formed near the outer, secondary ring at 0.326 Myr. The velocity field is plotted with respect to the local Keplerian velocity, and the vortex can be seen in the region marked with the black circle. Vortex formation can be triggered by Rossby wave instability, which can develop at pressure bumps near a sharp viscosity transition in a protoplanetary disk (Lovelace et al. 1999). Such vortices are capable of collecting dust and thus aid formation of planetesimals if they are long-lived. We found that the vortices formed near azimuthal asymmetries in the gas surface density in a ring and were usually in proximity with the gravitationally unstable regions occurring in the outermost parts of a ring structure. The vortices were always cyclonic (i.e., rotating in the same direction as the disk) and hence naturally unstable in a Keplerian disk owing to strong shear (e.g., Godon & Livio 1999). Regály & Vorobyov (2017) showed that the self-gravitating vortices are subject to strong torques and hence are prone to a significantly faster decay as compared to non-self-gravitating ones. We confirmed the rapid decay of the vortices by restarting parts of the simulation and obtaining a high temporal output of 20 yr. The vortices completely disappeared within this short time of disk evolution; however, note that their dissipation at 4 au should proceed quicker, over the local dynamical timescale. Their proximity to GI fragmentation suggests that the instability formed spiral arms and the resulting waves in turn generated local, short-lived vortices. One caveat of our simulations is that the vortices were relatively small and hence only marginally resolved—since the grid resolution at 4 au was about 0.08 au in the radial direction. However, the vortices were cyclonic and formed near locally unstable regions, and both of these factors contribute toward a rapid decay. Thus, additional grid resolution should not change the overall results of the simulations.

In order to study the rings formed in the layered disk model quantitatively, we define a ring to be the structure corresponding to the location of maximum surface density at a given time. From the spacetime diagrams in Figure 4 we can see that a ring typically started its evolution when there was another, more prominent ring closer to the star. This phase seemed to be about half of the total lifetime of a ring. Thus, with the above definition we cannot track the entire evolution of a ring; however, we should obtain estimates that are accurate up to a factor of two. In Figure 7 we consider the evolution of the fiducial layered disk model between 0.30 and 0.46 Myr when the rings are fairly stable in time. The first row shows the location of the ring, where the inward drift can be clearly noticed, as well as large discontinuities occurring at about 0.305 and 0.375 Myr. These discontinuities are associated with the instability of the most prominent ring and correspond to FUor-like outbursts. The vertical dashed lines mark these accretion events. Neglecting the first 0.2 Myr when the ring structures were highly asymmetric and transient, we calculated the average rate of migration of the rings to be about −25 au Myr⁻¹ (depicted by the diagonal red dashed line), while the average duration between the MRI outbursts was about 38,000 yr.

Zoom In Zoom Out Reset image size

The next two rows of Figure 7 show the surface density and α_eff at the location of the rings. As a ring formed a bottleneck for mass transport within the dead zone, the incoming material from the outer disk was accumulated in its vicinity. Thus, as the ring evolved, its surface density increased almost monotonically. The residual viscosity in the dead zone, which is induced by the action of the active surface layers, is parameterized via α_d (Equation (8)). Thus, as Σ_d increased, the α_eff decreased. The local value of α_eff could jump over two orders of magnitude during the MRI instability. Note that we did not catch the actual accretion event occurring at 0.305 Myr owing to the sparse rate of output (every 1000 yr), but we did for the outburst at 0.375 Myr. The fourth row of Figure 7 shows the evolution of the midplane temperature at the ring location. The temperature exceeded T_crit = 1300 K at the discontinuity at 0.375 Myr, indicating that the instability in the rings was indeed caused by MRI activation.

The last two rows of Figure 7 show the kinematic viscosity (ν_kin = $\alpha {{c}_{s}}^{2}$/ Ω_k) and the viscous timescale (T_visc = r²/ν_kin) at the ring location. As the ring evolved, the kinematic viscosity at its location decreased owing to the decrease in α_eff, indicating a narrower bottleneck for mass transport. As expected, the ν_kin briefly increased to a large value, and the bottleneck in mass transport disappeared at the accretion events. The viscous timescale of the rings stayed approximately constant at about 1 Myr. This implies that if left alone, the rings themselves would evolve and dissipate over this timescale. We hypothesize that the inward migration of a ring was primarily facilitated by large-scale spiral waves induced by GI, resulting in a rate of migration faster than the viscous timescale. The lifetime of a typical ring was much shorter than the viscous timescale because MRI was always triggered before the rings could show viscous dissipation. When MRI was triggered, the viscous timescale at the ring location decreased significantly, to about 1000 yr, for the outburst near 0.375 Myr. The material was accreted within this very short time as compared to the average viscous timescale of the rings. We found that the duration of the luminosity outburst (luminosity was sampled at a higher frequently) was in agreement with this minimum value of the viscous timescale. Notice that in Figure 4 as the simulation progressed the inward migration time of successive rings increased. We found that the asymmetries in the rings and GI fragmentation, both of which generated spiral waves, could be induced by strong perturbations in the inner disk—e.g., an infalling clump or an MRI outburst (described below). As the simulation evolved, the frequency of the clumps decreased in general (see Figure 4, Q-parameter), which should result in fewer GI spirals and the observed slower migration rate of the rings.

With the understanding of the properties and dynamical nature of the inner disk, we can propose a coherent picture of how the rings form and evolve in the layered disk models. As the disk evolves, the gas does not accumulate at the outer edge of the dead zone because it undergoes a gradual transition, as well as the action of transient GI spirals across the dead zone. When the gas is channeled in the inner disk, the viscous torques work in the direction of reinforcing the pileup of material. The positive feedback between viscous torques and accumulated gas produces the ring structures. Since the rings are marginally gravitationally unstable, perturbations of the inner disk may give rise to GI fragmentation and spiral waves. The associated transport of angular momentum results in the inward migration of the rings. The rings form within a few astronomical units from the inner edge of the dead zone; however, the exact location at a given time would depend on the interactions of all of the above factors. The MRI instability is ultimately triggered in the disk, and a ring rapidly depletes the accumulated material onto the star, resulting in an outburst.

The dynamical gaseous rings formed in a magnetically layered protoplanetary disk lie in the region that is of greatest interest for terrestrial planet formation. The dust grains that exist in the molecular cloud need to grow by about 13 orders of magnitude in order to form planetesimals (Williams & Cieza 2011). In this section we investigate the implications of such ring structures on the accumulation and growth of dust. In Figure 8 we compare the inner disk structure of the fiducial layered disk model with that of the fully MRI-active disk with respect to quantities that are relevant for the accumulation of dust—vertically integrated pressure, dust grain fragmentation radius (a_frag), and grain size at Stokes number unity (a_St = 1), as well as Stokes number at the dust grain fragmentation radius (St). The leftmost column in this figure shows the azimuthal profiles, while on the right we compare the spatial distributions of these quantities in the inner 10 × 10 au box. These quantities are calculated at 0.35 Myr, i.e., they correspond to the snapshot of the inner disk presented in Figure 3.

Zoom In Zoom Out Reset image size

Consider the first row of Figure 8, which compares the vertically integrated pressure for the two fiducial models. The azimuthally averaged pressure for the fully MRI-active model was monotonically decreasing with radial distance, although it showed a weak spiral structure within. We know that dust particles experience a drag due to their relative velocity with respect to the gas in the disk in the direction of the pressure gradient (Weidenschilling 1977). Thus, in the case of a fully active disk, the dust is predicted to drift toward the star. A weak spiral structure and associated local pressure maxima were shown to be inefficient in trapping dust particles in fully MRI-active disks, perhaps with the exception of the corotation region of the spiral pattern with the disk (Vorobyov et al. 2018). In contrast, the layered disk model showed several maxima in the azimuthally averaged pressure near the location of the gaseous rings. The dust grains can accumulate and grow in these regions of local pressure maxima.

The second row of Figure 8 compares the dust grain fragmentation radius for the two fiducial models, calculated as

$$\begin{eqnarray}&&{a}_{\mathrm{frag}}=\displaystyle \frac{2{\rm{\Sigma }}{v}_{\mathrm{frag}}^{2}}{3\pi {\rho }_{{\rm{s}}}{\alpha }_{\mathrm{eff}}{c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 16 }$$

where the assumed value of fragmentation velocity is v_frag = 10 m s⁻¹ and the internal density of the dust aggregate is ρ_s = 1.6 g cm⁻³. The a_frag is a measure of the maximum size that the dust particles can achieve before they are destroyed through the process of fragmentation, which tends to limit the dust size in the inner disk (Birnstiel et al. 2012; Vorobyov et al. 2018). The Stokes number for a particle of size a, under typical assumptions, is calculated as

$$\begin{eqnarray}&&\mathrm{St}=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{k}}}{\rho }_{{\rm{s}}}a}{\rho {c}_{{\rm{s}}}},\end{eqnarray} \tag{ 17 }$$

where the volume density $\rho ={\rm{\Sigma }}/\sqrt{2\pi }H$, with H being the local vertical scale height. The Stokes number describes the aerodynamic properties of dust grains and thus the coupling of the dust to the gas flow. The azimuthal plot of a_frag also shows the grain size at each radius that yields the Stokes number of unity. Dust particles at this size will achieve the maximum radial drift velocity and therefore show the quickest radial migration. Note that dust particles do not necessarily reach this size and barriers such as fragmentation, drift, bouncing, etc., tend to limit their growth, unless these are overcome by means of other mechanisms, e.g., streaming (Youdin & Goodman 2005) or gravitational (Boss 1997) instability. In the case of a fully MRI-active disk, the dust particles encountered a fragmentation barrier when they were a few millimeters across. However, in the layered disk model the fragmentation barrier was much larger in comparison because of the large local surface density and low α_eff. Within the dead zone the maximum particle size remained mostly between a few centimeters to 10 m, while occasionally approaching about 100 m in the vicinity of the rings.

The third row in Figure 8 shows the Stokes number of the dust particles at the dust grain fragmentation radius. The small values of Stokes number for the fully MRI-active disk mean that the dust particles were mostly coupled to and moved with the gas. Since the dust concentration efficiency is highest for a Stokes number of unity, the dust particles with this size (a_frag) would accumulate very quickly (over the local orbital timescale) in the vicinity of the rings. Numerical studies of gas and dust evolution in a fully MRI-active protoplanetary disk show that the radial inward drift of dust particles poses a significant barrier for planet formation (Birnstiel et al. 2010; Vorobyov et al. 2018). The gaseous rings observed in our simulations could naturally form dust traps, effectively preventing this radial drift and providing sites for the accumulation of dust particles. However, the dynamic nature of the rings complicates the situation for their applicability for planet formation. As we saw in Section 3.2, although the viscous timescale inside the rings was sufficiently long (≈1 Myr), the lifetime of the rings was much shorter—only a few percent of this value—before the MRI instability was triggered. Evolution of the dusty component needs to be taken into account in order to estimate the timescale of the dust growth in this atypical environment and to determine whether the rings can successfully halt the radial drift long enough to form planetesimals.

So far we concentrated on the differences between a typical layered disk model and the classical approach of a fully MRI-active protoplanetary disk. As discussed in Section 2.3, the parameters used within our layered disk model are constrained but not certain. The uncertainty in T_crit arises owing to the complexity of the dust and gas physics, as well as that in estimating when the disk can be considered MRI-active. On the other hand, Σ_a may depend on the local conditions in the disk such as the flux of ionizing radiation, strength of magnetic field, and properties of dust grains. In this section we investigate the dependence of the disk structure and evolution on these model parameters. We present model1_T1500_S100 and model1_T1300_S10 to demonstrate the effects of changing critical temperature (T_crit = 1500 K) and active layer thickness (Σ_a = 10 g cm⁻²), respectively. Most protostars have a mass less than 1 M_⊙, so the dependence on mass is demonstrated with the help of model2_T1300_S100. In this case the initial cloud core mass was 0.346 M_⊙, which is about 30% of the mass of the core in the fiducial model. As described in Table 1, these models are otherwise identical to the fiducial layered disk model. In the interest of conciseness, the above three models are presented in detail, while the remaining simulations are referred to only when required.

Figure 9 shows the typical disk structure in terms of azimuthal profiles of the three models—model1_T1500_S100, model1_T1300_S10, and model2_T1500_S100. We considered only Σ and α_eff, as these two quantities are sufficient to demonstrate the salient properties of the rings and the dead zone. Note that the abscissa shows the inner 200 au region in logarithmic units. These stationary profiles are obtained at 0.35 Myr for the higher-mass models (model1) and at 0.18 Myr for the lower-mass model (model2). Here we describe the main features of the typical ring structures, while the underlying causes between the differences are explained later in this section. Consider model1_T1500_S100 (red lines) in Figure 9 showing the effects of increasing T_crit, from 1300 to 1500 K on the ring structures. Three clearly formed rings can be seen in this case as compared to only two for the fiducial model at this time. The maximum surface density and minimum in α_eff at the location of the most prominent ring were similar to those of the fiducial model. The outer extent of the dead zone—as shown by the red dashed line in the bottom panel—was about 16 au, almost equal to that of the fiducial model. The effect of changing Σ_a from 100 to 10 g cm⁻² (model1_T1300_S10; blue lines) was much more significant. Although the maximum value of Σ was not affected, a much deeper dead zone in terms of α_eff was formed. Unlike the other models with Σ_a = 100 g cm⁻², the individual rings in this case were not as prominent and clearly distinguishable from the accumulated gas in the dead zone. The variations between the maximum and the minimum at a given radius showed larger perturbations, indicating relatively more azimuthal asymmetry. Also note that several maxima in gas surface density and consequently in pressure were formed, which could trap dust at several radii throughout the dead zone. The extent of the dead zone significantly increased to about 78 au, which is almost 5 times larger as compared to the fiducial model. Consider model2_T1500_S100 (green lines) depicting the effect of lowering the mass of the parent core. Only a single ring was formed at this time, while the maximum in Σ or the minimum in α_eff at the location of the ring were not significantly different as compared to the fiducial model. The extent of the dead zone was much smaller, only 3.3 au, which is about 5 times smaller than the fiducial model. Note that for this lower-mass model the relatively flat surface density reached Σ_a = 100 g cm⁻² at a much larger radius.

Zoom In Zoom Out Reset image size

In order to demonstrate the effects of model parameters on the temporal evolution, the spacetime diagrams of the three simulations are shown in Figures 10–12. These figures are similar to Figure 4 for the fiducial models. Consider Figure 10 for model1_T1500_S100, showing the effects of increasing T_crit over about 0.5 Myr. Notice that the qualitative behavior of the gaseous rings appears similar to the fiducial layered disk model. The rings originated at the distance of a few astronomical units and migrated inward, which resulted in the diagonal patterns in the gas surface density. The rate of inward migration was similar to the fiducial model. The rings also disappeared suddenly, which indicates an instability similar to that described earlier, resulting in FUor-like accretion events. The overall extent of the dead zone, denoted by the black contour in α_eff, was also similar to the fiducial model throughout the evolution. However, compared to the fiducial model, the ring structures display more complex behavior that was difficult to quantify. As shown earlier, more rings were formed simultaneously, e.g., between about 0.35 and 0.40 Myr, there are three rings in the disk as compared to two for the fiducial model. Another example of complexity is that the two rings merged to form a single ring at about 0.25 Myr. With the increase in T_crit, the MRI would be triggered at a higher temperature and the gas in the disk would remain in the disk for a longer time. Since T_crit was higher, the rings could reach a smaller radius before MRI instability was triggered. Thus, the inner rings were long-lived, and gas moving in from the outer disk formed multiple outer rings. The average time between MRI eruptions was about 55,000 yr, which is longer than the fiducial model and thus consistent with late onset of the instability. The complexity of the rings can also be seen in the spacetime diagram of α_eff. Thus, we conclude that an increased T_crit rendered the inner disk structure more complex in terms of the ring structures; however, the overall behavior was not affected significantly. We confirmed this conclusion with the comparison of the remaining three simulations with T_crit =1500 K against their counterpart with T_crit = 1300 K (see Table 1). This result can be expected since T_crit is well constrained and was increased by only about 15%.

Zoom In Zoom Out Reset image size

The spacetime diagrams for model1_T1300_S10 are presented in Figure 11, showing the effects of a decrease in the active layer thickness. Consider the behavior of gas surface density in the top panel. Similar to the fiducial model, the ring structures are formed in the inner disk at a radius of a few astronomical units. As shown in Figure 9, these rings did not display a significant enhancement of gas surface density in the rings as compared to the surrounding material. The major difference in the temporal behavior—as compared to the previous simulations—was that the rings showed ∼0.02 Myr long diagonal patterns, and then no significant evolution. Each of the inward migrations was preceded by an MRI accretion event that was qualitatively similar to those described in Section 3.2, and can be seen as the discontinuities in Σ, as well as jumps in α_eff, in Figure 11. Since the rings in model1_T1300_S10 were not as pronounced as in the previous cases, instead of complete disappearance of a ring, a partial "cave-in" of the inner disk region occurred. During an MRI outburst, the inner disk was rapidly depleted of material, and the newly formed ring migrated inward owing to the action of GI spirals. Notice that α_eff fell below 10⁻⁵ in the vicinity of the rings, which was an order of magnitude lower as compared to the fiducial model. We conjecture that once the rings migrated inward they were stable against fragmentation by inefficient cooling, and their slow evolution onward was a direct consequence of deepening of the rings with respect to α_eff. From Equation (8) it can be deduced that if Σ_a is decreased, α_rd will also decrease accordingly. Thus, lowering Σ_a by an order of magnitude resulted in a proportional decrease in α_eff in the region of the rings. The black contour in the bottom panel shows a much larger extent of the dead zone (as compared to the simulations with Σ_a = 100 g cm⁻²), about 78 au, which lasted throughout the 0.5 Myr of evolution. The width of the low-viscosity region (α_eff = 10⁻⁴) was also increased with respect to the fiducial model because of this accumulation of gas. As a consequence, the viscous timescale in the vicinity of the rings was of the order of 10 Myr, which is again 10 times larger than the fiducial model. Due to the much wider dead zone and lower viscosity, the region must form a much more efficient bottleneck for mass and angular momentum transport as compared to the fiducial model. The higher surface gas density in the dead zone region possibly translated into strong GI spirals and a faster initial migration of the rings. These spirals may have resulted in a larger scatter between the maxima and minima seen in Figure 9. The other three models with Σ_a = 10 g cm⁻² from Table 1 evolved in a similar manner, showing a significantly slower evolution of the rings, as well as a much deeper and wider region with respect to α_eff.

Zoom In Zoom Out Reset image size

Figure 12 depicts the evolution of model2_T1300_S100, where we show the effects of a decrease in the mass of the cloud core. The change in the mass available for disk formation significantly affected the evolution of the inner disk region, which includes the gaseous rings as well as the dead zone. Notice that with a smaller cloud core mass, the disk started forming at an earlier time at about 0.025 Myr, as compared to about 0.08 Myr for the higher-mass models. The spacetime diagram for Σ shows that a ring started to form after some delay at about 0.15 Myr. As expected, the ring migrated inward; however, this "ring phase" (approximately defined by the span of α_eff = 10⁻⁴ contour) did not last long, and the ring disappeared within only about 0.05 Myr. As compared to the fiducial model, this region of low α_eff was much smaller in width as well. We hypothesize that initially the ring formed because the gas surface density could be maintained above Σ_a owing to the material moving in from the collapsing core through the outer disk. Due to the cloud core being low mass, the outer disk contained a relatively limited amount of mass. As the gas was transported inward through the active layer, it was not replenished from the outer disk. Thus, a layered structure could not be maintained for an extended period of time. Despite the low mass of the disk, we observed GI fragmentation and spirals in the rings. This may be because the low mass of the star increased dynamical time and the increased cooling efficiency was sufficient for the onset of GI. The dead zone, as depicted by the black contour, also disappeared within a short time at about 0.25 Myr, and the disk essentially became fully MRI-active thereafter. Note that the yellow contour in the top panel showing Σ = Σ_a grossly overestimated the extent of the dead zone in both radius and time. There were no MRI-triggered instabilities throughout the evolution, including when the ring was fully terminated. This indicates that the midplane temperature never reached T_crit in the relevant region and the ring eventually viscously dissipated. Thus, we conclude that the mass of the initial cloud core can significantly affect the structure and evolution of the inner disk.

Zoom In Zoom Out Reset image size

In order to show the effect of Σ_a on the entire evolution of the ring phase, we investigate model2_T1300_S10. Figure 13 depicts the spacetime diagrams for this model, showing the effect of a low active layer thickness on an initially low-mass cloud core. As expected, the dead zone formed was much wider and the rings formed were deeper in terms of α_eff, as compared to the analogous model with a higher Σ_a (model2_T1300_S100). However, the duration of the ring phase was several times longer in comparison. This model also showed an absence of MRI instabilities and associated outbursts, although the dead zone was maintained throughout the evolution. Thus, we conclude that the inner disk structure and the properties of the rings were highly dependent on Σ_a and the initial cloud core mass but only marginally dependent on T_crit.

Zoom In Zoom Out Reset image size