TWO-DIMENSIONAL SIMULATIONS OF FU ORIONIS DISK OUTBURSTS

The FU Orionis systems are a small but remarkable class of young stellar objects, which undergo outbursts in optical light of 5 mag or more (Herbig 1977). While the rise times for outbursts are usually very short (∼1 − 10yr), the decay timescales range from decades to centuries. In outbursts, they exhibit F-G supergiant optical spectra and K-M supergiant near-infrared spectra dominated by deep CO overtone absorption. The FU Ori objects also show distinctive reflection nebulae, large infrared excesses of radiation, wavelength-dependent spectral types, and "double-peaked" absorption line profiles (Hartmann & Kenyon 1996). The accretion disk model for FU Ori objects proposed by Hartmann & Kenyon (1985, 1987a, 1987b) and Kenyon et al. (1988) can explain the peculiarities enumerated above in a straightforward manner.

Bell & Lin (1994) suggested that the thermal instability (TI) was responsible for FU Ori outbursts. To agree with the outburst mass accretion rates (∼10⁻⁴ M_☉ yr⁻¹) and the outburst duration time (∼100 yr), with the TI, Bell & Lin used very low viscosities both before and after the outbursts. However, this theory predicts that the outburst is confined to inner disk radii ≲0.1 AU, whereas Zhu et al. (2007, 2008) found that the high accretion rate region must extend to ≳0.5 AU, based on modeling Spitzer IRS data.

Bonnell & Bastien (1992) suggested that a binary companion on an eccentric orbit could perturb the disk and cause repeated outbursts of accretion. However, in the case of FU Ori, no close companion is known, and no radial velocity variation has been detected (Petrov & Herbig 2008). Vorobyov & Basu (2005, 2006) suggested that gravitational instabilities might produce the non-steady accretion, but could not carry out the calculation within the inner disk; thus, the details of the accretion event are uncertain. Lodato & Clarke (2004) suggested that planets might dam up disk material until a certain point, resulting in outbursts. This model also depends on the TI and assumes a very low viscosity as in Bell & Lin (1994).

The rediscovery and application of the magnetorotational instability (MRI) to disk angular momentum transport have revolutionized our understanding of accretion in ionized disks (e.g., Balbus & Hawley 1998 and references therein). However, it was recognized that T Tauri disks are cold and thus likely to have such low ionization levels that the MRI cannot operate, at least in some regions of the disk (e.g., Reyes-Ruiz & Stepinski 1995). Gammie (1996) suggested that non-thermal ionization by cosmic rays could lead to accretion in upper disk layers, with a "dead zone" in the midplane. This picture implies that accretion is not steady throughout the disk, and the mass should pile up at small radii.

Armitage et al. (2001) considered the time evolution of layered disks including angular momentum transport by the gravitational instability (GI). Within a certain range of mass addition by infall in the outer disk, they obtained episodic outbursts of accretion at ∼10⁻⁵ M_☉ yr⁻¹ lasting ∼10⁴ yr. However, these outbursts are too weak and last too long compared with FU Ori outbursts. Gammie (1999) and Book & Hartmann (2005) were able to get outbursts more like FU Ori behavior with a different set of parameters. The Armitage, Gammie, Book, & Hartmann models were vertically averaged and had schematic treatments of radial energy transport.

Earlier models of layered disk evolution have been based on phenomenological one-zone, or two-zone, models for the vertical structure of the disk. It is quite difficult to formulate these models in a consistent way, particularly near transition radii where the temperature or ionization structure of the disk changes sharply and radial energy transport may be important. To make progress, it is necessary to consider more accurate vertical structure and more accurate energy transport models.

In this paper, we construct two-dimensional (axisymmetric) layered disk models. We find that the MRI triggering by the GI can explain FU Ori-type outbursts if the MRI is thermally activated at T_M ∼ 1200 K. In Section 2, the equations and methods of the two-dimensional models are introduced. In Section 3, outburst behavior is described. In Section 4, we focus on vertical structure and convection in the disk. We discuss implications of the two-dimensional simulations compared with observations in Section 5, and summarize our results in Section 6.

In this work, we are particularly interested in accurately modeling the flow of energy within the disk. A full treatment would numerically follow the development of MRI, GI, and radiative, convective, and chemical energy transport in three dimensions. As this is not yet possible, we use the phenomenological α-based viscosity to represent the angular momentum transport by the MRI and GI and heating by dissipation of turbulence. An accurate treatment of energy balance (viscous heating and radiative cooling) is essential to understand the outburst, especially for the TI. We implemented the realistic radiative cooling in the ZEUS-type code of McKinney & Gammie (2002).

To test our code, we set a hot regioin an optically thick disk region and then let the radiative energy diffuse without evolving the hydrodynamics. The hot region expands nicely in a round shape in R–Z space with the correct flux. We have also tested the total energy conservation during our outburst simulations, finding that the energy deficit during outburst is only <1% of the total energy radiated.

2.1. Viscous Heating

The momentum equation and energy equations are

$$\begin{equation} \rho \frac{D v}{Dt}=-\nabla p-\rho \nabla \Psi -\nabla \cdot \Pi \,, \end{equation} \tag{ 1 }$$

$$\begin{equation} \rho \frac{D}{Dt}\left(\frac{e}{\rho }\right)=-p \nabla \cdot \nu +\Phi -\Lambda \,. \end{equation} \tag{ 2 }$$

Here, as usual, D/Dt ≡ ∂/∂t + ν · ∇ is the Lagrangian time derivative, ρ is the mass density, ν is the velocity, e is the internal energy density of the gas, p is the gas pressure, Ψ is the gravitational potential, and Λ is the radiative cooling rate. The dissipation function Φ is given by the product of the anomalous stress tensor Π with the rate-of-strain tensor ξ,

$$\begin{equation} \Phi =\Pi _{ij}\xi ^{ij} \end{equation} \tag{ 3 }$$

(sum over indices), where the anomalous stress tensor is the term-by-term product

$$\begin{equation} \Pi _{ij}=2\rho \nu _{V} \xi _{ij}S_{ij} \end{equation} \tag{ 4 }$$

(no sum over indices), where S_ij is a symmetric matrix filled with 0 or 1 that serves as a switch for each component of the anomalous stress. In current simulations, we switched on all the component of the anomalous stress. It is not yet known whether any prescription for Π_ij is a good approximation for modeling magnetohydrodynamical (MHD) turbulence. We set the viscosity coefficient to ν_V = α(c²_s/Ω)sin^3/2θ, which vanishes rapidly near the poles (McKinney & Gammie 2002). α is the viscosity parameter and Ω is the angular velocity at radius R.

The equation of state is

$$\begin{equation} p=\left(\gamma -1\right)e\,, \end{equation} \tag{ 5 }$$

and the sound speed is

$$\begin{equation} c_{s}^{2}=\gamma \frac{k_{B}T}{\mu }, \end{equation} \tag{ 6 }$$

where γ and μ are the adiabatic index and the mean molecular weight. We assume γ = 5/3 and μ = 2.4M_P, where M_P is the mass of one proton. The adiabatic index and mean molecular weight will change at different temperatures and pressure (especially during outburst when the disk is fully ionized) and a fully self-consistent treatment would use a tabulated equation of state. We will discuss a more realistic equation of state in a later paper.

2.2. Radiative Cooling

In an accretion disk, radiation is the main cooling mechanism, and radiation energy and pressure can affect the disk dynamics. Numerically, there are two different approaches to treat the radiative transport in an accretion disk. One is to solve the full radiative transfer equation

$$\begin{equation} \rho \frac{D}{Dt}\left(\frac{E}{\rho }\right)=-\nabla \cdot F-\nabla \nu:P +4\pi \kappa _{P}B-c\kappa _{E}E \end{equation} \tag{ 7 }$$

with the material energy equation

$$\begin{equation} \rho \frac{D}{Dt}\left(\frac{e}{\rho }\right)=-p\nabla \cdot \nu-4\pi \kappa _{P}B+c\kappa _{E}E \end{equation} \tag{ 8 }$$

in both the optically thick and thin regions by an implicit method (Turner & Stone 2001; Wünsch et al. 2005; Hirose et al. 2006). Here, as above, ρ, e, ν, and p are the gas mass density, energy density, velocity, and scalar isotropic pressure, respectively, while E, F, and P are the total radiation energy density, momentum density, and pressure tensor, respectively. κ_P and κ_E are Planck mean and energy mean absorption coefficients. An implicit method is required because of the extremely short time step needed by the explicit method (grid size over the speed of light), but the implicit method is time consuming, requiring solution of an N² × N² matrix equation at each time step.

A second, less accurate but less costly, approach is to add a radiative source (heating/cooling) term in the energy equation, use the diffusion approximation in the optically thick region, and use an approximation to treat the optically thin region (Boss 2002; Boley et al. 2006; Cai et al. 2008). When the radiation pressure is negligible, as in the case for protostellar disks, and the radiation field is nearly in equilibrium (D(E/ρ)/Dt = 0), Equations (7) and (8) become

$$\begin{equation} \rho \frac{D}{Dt}\left(\frac{e}{\rho }\right)=-p\nabla \cdot \nu-\Lambda, \end{equation} \tag{ 9 }$$

where Λ = ∇ · F and F is the radiative flux.

We examine the justification for ignoring the advection term (within the term on the left side of Equation (7)) and the term describing work done by radiation on the gas (second term of the right-hand side of Equation (7)) following the criterion given by Mihalas & Weibel Mihalas (1984). If in one system, L is the characteristic size, u is the characteristic velocity, and τ is the optical depth of L, the terms describing the advection and the work done by the radiation can be ignored if τ≫1 and τ × u/c≪1, which is called static diffusion limit. This criterion can be simply understood as follows.

The timescale for radiation to propagate across a region of length L and mean optical depth τ is of order

$$\begin{equation} t_E \sim L \tau / c \,. \end{equation} \tag{ 10 }$$

For comparison, the material crossing time of this region is

$$\begin{equation} t_u \sim L / u \end{equation} \tag{ 11 }$$

so that the ratio of these two timescales is

$$\begin{equation} {t_E \over t_u} \sim {u \over c} \tau \,. \end{equation} \tag{ 12 }$$

If this ratio is far less than 1, the region is in radiative quasi-equilibrium before the material carrying radiation energy away so that the energy emitted and absorbed equals to the diffusion term (Equation (7)). In a three-dimensional accretion disk, this typical velocity is the Keplerian velocity. However, in our two-dimensional, axisymmetric simulation, this velocity is just the mean radial velocity of the gas. Thus, t_u is close to the viscous time which is longer than the time for sound to cross a scale height H by a factor of order α⁻¹(R/H)².

In our simulations, even at the maximum temperatures achieved in the innermost disk, u/c ≲ 10⁻⁵ (Section 4) and so the approximation of radiative equilibrium should be valid except at very large optical depths. It turns out that, during the outburst, vertical Rosseland mean optical depths ∼10⁴ in the innermost disk. Thus, our approximation is justified, though it may break down at smaller radii than we have considered here.

In detail, the radiative transfer in the optically thick disk, defined as the region where τ>τ_thick, is given in the diffusion approximation by

$$\begin{equation} \Lambda =\triangledown \cdot F, \end{equation} \tag{ 13 }$$

where the flux F is

$$\begin{equation} F=-\frac{16\sigma \beta T^{3}}{\chi _{R}\rho }\triangledown T, \end{equation} \tag{ 14 }$$

and β is the flux limiter,

$$\begin{equation} \beta =\frac{2+y}{6+3 y+y^{2}}\,, \end{equation} \tag{ 15 }$$

with

$$\begin{equation} y=\frac{4}{\chi _{R}\rho T}\mid \!\triangledown T \!\mid \,, \end{equation} \tag{ 16 }$$

where χ_R is the Rosseland mean of the absorption coefficient.

The optically thin region is more complicated. In this region, we restrict our considerations to radiative transport in the vertical direction only. We have developed two methods to solve this problem in collaboration with members of the Indiana disk group (D. Durisenet al. 2008, private communication). The first method considers the cooling rate of every cell in the optically thin region can be written as

$$\begin{equation} \Lambda =\rho \kappa _{P}(T)(4\sigma T^{4}-4\pi J(Z))-\rho \kappa _{P}(T_{e})\sigma T_{e}^{4}{\rm e}^{-\tau _{e}(Z)}\,, \end{equation} \tag{ 17 }$$

where κ_P(T) is the Planck mean opacity at the temperature T, T_e is the temperature of the envelope, and also J(Z) and τ_e(Z) are the mean intensity and the optical depth of the envelop at height Z in the disk atmosphere. The first term is the energy emitted by the cell, while the second and third terms are the energy absorbed from the optically thick region and from the envelope. Using the two-stream approximation,

$$\begin{equation} J(Z)=\frac{1}{4\pi }\int I {\rm d}\Omega =\frac{1}{4\pi }\left(2\pi I_{+}+2\pi I_{-}\right). \end{equation} \tag{ 18 }$$

This treatment assumes flux conservation and thus ignores energy dissipation in the atmosphere. The net flux F_b from the disk photosphere is then

$$\begin{equation} F_{b}=\int I \, \mu \, {\rm d}\Omega =\pi \left(I_{+}-I_{-}\right)\,, \end{equation} \tag{ 19 }$$

where μ is the cosine of the angle measured from the vertical. Thus,

$$\begin{equation} J(Z)=\frac{1}{2}\left(I_{+}+I_{-}\right)=\frac{F_{b}}{2\pi }+I_{-}(Z). \end{equation} \tag{ 20 }$$

F_b is approximately

$$\begin{equation} F_{b}=\frac{4\sigma \left(T_{b}^{4}-T_{e}^{4}\right)}{3\left(\tau _{b}+2/3\right)}, \end{equation} \tag{ 21 }$$

where T_e is a term which accounts for the external irradiation (see Equation (10) in Boley et al. (2006)). T_b and τ_b are the temperature and optical depth of the first optically thick cell right below the photosphere (again, see Boley et al. (2006)). Finally, we must set I₋(Z) to obtain a closed expression for Λ:

$$\begin{equation} I_{-}(Z)= \int _{Z}^{\infty }\rho \kappa _{P}\frac{\sigma T(z)^{4}}{\pi }{\it d}z+\frac{\sigma T^{4}_{{\rm env}}{\it e}^{-\tau _{en}(Z)}}{\pi }\,, \end{equation} \tag{ 22 }$$

which is the intensity from the optically thin region above height Z. Theoretically, the first term should also have e^−τ. However, since we only consider the optically thin region here, we ignore it to simplify the calculation and save computational time.

The second method is much simpler but limited to this particular problem with no irradiation. As we do not consider irradiation, we can use the usual gray atmosphere result. Without the third term in Equation (17), we have

$$\begin{equation} \Lambda =\rho \kappa _{P}(T)(4\sigma T^{4}-4\pi J(Z))\,, \end{equation} \tag{ 23 }$$

and in a gray case

$$\begin{equation} J=\frac{\sigma }{\pi } T^{4}(Z)=\frac{\sigma }{\pi }\left(\frac{3}{4}T_{{\rm eff}}^{4}\left(\tau +\frac{2}{3}\right) \right)\,, \end{equation} \tag{ 24 }$$

where T_eff is given by the boundary flux F_b = σT⁴_eff as Equation (21). In detail, counting θ grid cell index i as increasing with increasing θ and τ, we set the first grid with τ(R, θ_i)>3 as the boundary grid and τ_b = τ(R, θ_i). Then, we derive the effective temperature (or equivalently, the boundary flux) using the gray approximation

$$\begin{equation} F_b=\sigma T_{{\rm eff}}^4 = \frac{4\sigma T_{b}^{4}}{3(\tau _{b}+2/3)} \,. \end{equation} \tag{ 25 }$$

With T_eff, we use Equation (24) to set the temperatures of all cells at R, θ < θ_i, and use the diffusion Equation (14) at all r, θ ⩾ θ_i. This way of matching the "optically thin" (gray) atmosphere to the "optically thick" (diffusion equation) interior removes the temperature jump seen at the boundary between the two treatments seen in the simulations of Cai et al. (2008).

Both methods give similar results if there is no irradiation (see Appendix A)—basically, they are the same except for the assumed directionality of the radiation field. We use the second method in this paper, although the first method is of interest for cases in which external irradiation is important.

We have adopted the Bell & Lin (1994) opacity fits to facilitate comparison with previous calculations (Bell & Lin 1994; Armitage et al. 2001); subsequent improvements in low-temperature molecular opacities will modify results at low temperatures (Zhu et al. 2009), but should not substantially affect the outburst behavior driven at higher temperatures.

2.3. Angular Momentum Transport

To complete the model, we need a prescription for the dimensionless viscosity α. The key physics we want to capture is the dependence of the MRI and GI transport on position in the disk (see, e.g., Miller & Stone 2000 for MRI activity as a function of height in the disk).

For the MRI, we set α_M = 0.1 when T > T_M, where T_M is a critical temperature (typically 1200 K), where thermal ionization makes the disk conductive enough to sustain the MRI.⁵ In addition, we calculate Σ₊ = ∫^θ₀Rdθ'ρ (and in models where the disk is not assumed to be symmetric about the midplane, Σ₋ = ∫^π_θRdθ'ρ) and set α_M = 0.1 if Σ₊ or Σ₋ < Σ_crit, where Σ_crit is the maximum surface density of the MRI active layer on one side of the disk (in our current model, we assume Σ_crit=50 g cm⁻²). Otherwise the disk is magnetically "dead" and α_M = 0. Also to avoid the strong MRI viscous heating at the disk's upper photosphere with very low density, we set α = 0 for the region with τ < 10⁻³.

The value of α_M predicted by theory is uncertain, with widely varying estimates depending on mean magnetic fluxes, magnetic Prandtl numbers, and other assumptions (Guan & Gammie 2008; Johansen & Levin 2008; Fromang et al. 2007). Given these uncertainties, observational constraints take on special significance. A recent review of the observational evidence by King et al. (2007) argues that α_M must be large, of the order 0.1–0.4, based in part on observations of dwarf novae and X-ray binaries where there is no question of GI. Our empirical study of FU Ori, where we were able to constrain the radial extent of the high accretion state, implies values of α ∼ 0.1. (Zhu et al. 2007).

For the GI, we calculate the Toomre Q parameter

$$\begin{equation} Q = {c_s \Omega \over \pi G \Sigma } \,. \end{equation} \tag{ 26 }$$

When Q ∼ 1 the disk is gravitationally unstable. Full global models (Boley et al. 2006) and physical arguments (Gammie 2001) suggest that if the disk is sufficiently thin angular momentum transport by GI is localized and can be crudely described by an equivalent α. We set

$$\begin{equation} \alpha _{Q}(R,\theta)={e}^{-Q(R)^{2}}, \end{equation} \tag{ 27 }$$

which is independent of θ. This has the sensible effect of turning on GI-driven angular momentum transport when Q < 2 and turning off GI-driven angular momentum transport when Q is large. In regions with both GI and MRI, we set α = α_M (in these regions α_M ≫ α_Q).

The maximum α_Q is <0.1 during the simulations, consistent with our assumption that fragmentation does not occur. Our prescription for α_Q on Q is heuristic and other forms are plausible. Gammie (2001) and Lodato & Rice (2004, 2005) found that the amplitude of GI depends on the cooling rate in addition to Q, and Balbus & Papaloizou (1999) argued that the angular momentum flux carried by GI is a global effect and the local viscosity description can only be valid near the corotation radius. However, using self-consistent cooling, Boley et al. (2006) showed that the α_GI agrees with the Gammie description very well, which assumes that the GI viscous heating is locally balanced by the cooling. Recent smooth particle hydrodynamics (SPH) simulations by Cossins et al. (2009) also showed that GI disks with tightly wound spiral arms can be approximated by a local treatment.

We use spherical polar coordinates and assume axisymmetry in all our simulations. To sample both the inner disk and the outer disk, we use a logarithmic radial grid (radial grids are uniform in log(R − R_in)) from 0.06–40 AU and a uniform θ grid in (0.5 ± 0.15)π. This gives adequate radial resolution for both the inner and outer disk, and spans a large enough range in θ to include the optically thick region of the disk. Parameters for all the two-dimensional simulations are listed in Table 1.

3.1. Initial and Boundary Conditions

All our disk models begin with the initial temperature distribution

$$\begin{equation} T(R,\theta) = 300 (R/{\rm AU})^{-1/2} \ {\rm K}, \end{equation} \tag{ 28 }$$

and ρ = Σ(R)/(H2π^1/2)exp(−R²θ²/(2H²)), H ≡ c_s/Ω from 0.4–100 AU⁶, so the disk is approximately in dynamical equilibrium. We set Σ(R) using a $\dot{M}=10^{-5}\, M_{\odot }\ {\rm yr}^{-1}$ constant α (α = 0.1) one-zone accretion disk model (see Equation (3) in Whitney et al. 2003). Although this is not an exact steady-state model, the disk relaxes rapidly at the beginning of the simulation.

We apply outflow boundary conditions (Stone & Norman 1992) to all boundaries. In some simulations, the disk is assumed to be symmetric about the midplane to reduce computational cost; for these half-disk simulations, we use a reflecting boundary condition (Stone & Norman 1992) at the midplane. The inner boundary is set to either 0.06, 0.1, or 0.3 AU to examine how different inner radii affect disk evolution. These parameters are listed in Table 1.

We follow the evolution of the disk for 10⁵ yr. Typically several outbursts occur during the simulation. After the first outburst, the subsequent outbursts are similar and do not depend on the initial condition. Thus, we use the disk density and temperature distribution after the first two outbursts as the new initial condition to study the third outburst for various one-dimensional and two-dimensional models. The disk mass for this initial condition is ∼0.5 M_☉.

3.2. Outbursts

In the protostellar phase, the disk is unlikely to transport mass steadily from ∼100 AU all the way to the star at an accretion rate matching the mass infall rate 10⁻⁶ to 10⁻⁴ M_☉ yr⁻¹ from the envelope to the outer disk. The outburst sequence we find here is qualitatively similar to that found by Armitage et al. (2001), Gammie (1999), and Book & Hartmann (2005). Mass added to the outer disk moves inward due to GI, but piles up in the inner disk as GI becomes less effective at small R. Eventually, the large Σ and energy dissipation leads to enough thermal ionization to trigger the MRI. Then the inner disk accretes at a much higher mass accretion rate, which resembles FU Orionis-type outbursts. After the inner disk has drained out and becomes too cold to sustain the MRI, the disk returns to the low state. During the low-state mass continuously accretes from large radius (or from an infalling envelope) causing matter to pile up at the radius where GI becomes ineffective yet again, leading to another outburst.

Figures 1–5 show the disk temperature and density structure at four consecutive stages of the outburst for our fiducial model. Figure 1 is shown in R–Z coordinate while Figures 2 and 3 are shown in log₁₀R–θ coordinates (R from 0.06 to 40 AU, θ from 1.1 to π/2). Since the GI, MRI and TI are taking place at quite different radii, most of our plots are in log₁₀R–θ coordinates so that all of them can be seen in one figure, although it is less intuitive than R–Z coordinates. The black solid curves plotted on top of the color contours delimit regions of differing α. Regions above the top solid curve limit the low-ρ photosphere with τ < 10⁻³ to α = 0. Between the top and next lower solid curve lies the MRI-active layer with Σ₊ < 50 gcm⁻¹ and α = α_M = 0.1. Regions below the two solid curves constitute the "dead zone" with only α_Q if any. The vertical-dashed lines at large radii represent α_Q = 0.025.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Panel (a) of Figures 1–3 shows the disk temperature structure before the outburst. The disk is too cool to sustain any MRI except in the active layer. The contours in Figures 1–3 show the outer disk (R > 1 AU) is geometrically thick, with a ratio of vertical scale height H to radius H/R ∼ 0.2, and massive, while the inner disk is thin with a low surface density. The inner disk is stable against GI and cannot transport mass effectively. The mass coming from the GI driven outer disk piles up at several AU. The GI thus gradually becomes stronger at several AU and the disk heats up.

The disk midplane temperature at 2 AU eventually reaches the MRI trigger temperature (T_M∼1200 K) and MRI driven angular momentum transport turns on. Panel (b) of Figures 1–3 shows the stage when the thermal MRI is triggered at the midplane around 2 AU. The thermal MRI active region are outlined by the semicircular α contour around 2 AU in these figures. After the MRI is activated, α increases from α_Q(2 AU) ∼ 0.01 to α_M = 0.1, which leads to a higher mass accretion rate. Afterward, the MRI active region expands both vertically and radially until the whole inner disk becomes fully MRI active.

As more and more mass is transferred to the inner disk, the temperature continues to rise. At T ∼ 4000 K the disk becomes thermally unstable and T rises rapidly to ∼2 × 10⁴ K. Panel (c) of Figures 1–3 shows the moment the TI is initiated in the innermost disk. With the temperature leap, the pressure of the TI active region also increases significantly. Thus, the TI active region expands and compresses the TI non-active region, leading to a lower density of the TI active region than non-active region (lower density at midplane around 0.1 AU in panel (c) of Figure 3).

Panel (d) of Figures 1–3 shows that after the TI is triggered in the innermost disk, the TI front travels outward to ∼0.3 AU until the entire inner disk is in a "high state." The hot inner disk (R < 0.2 AU) is puffed up to H/R ∼ 0.15. The mass accretion rate is 2 × 10⁻⁴ M_☉ yr⁻¹ at this stage. The central temperature has a sharp drop at 0.3 AU, beyond which the disk is in a thermally stable state at T ≃ 3000 K. Details of the TI are discussed in Appendix B. At this transition from the high to low state around 0.2 AU, some temperature and density fluctuations are present, which results from convection. We will discuss convection in detail in Section 5.

The TI high state (T_c > 10, 000 K, panel (d) of Figures 1–3) lasts for ∼10² yr, draining the inner disk (R < 0.5 AU) until it is not massive and hot enough to sustain the TI high state; it then returns to the TI low state. The inner disk may temporarily return to the high state due to MRI-driven mass inflow but eventually mass is drained from the MRI active region, it becomes MRI inactive, and the disk returns to the state shown in panel (a).

Figure 4 shows the vertically integrated surface density at the above stages. Before the outburst, the surface density peaks around 2 AU (long-dashed curve); this radius divides the inner, GI stable disk from the outer GI unstable disk. Once the MRI is triggered at 2 AU, it transports mass into the inner disk effectively. At the moment the TI is triggered (dotted curve), the surface density of the inner disk is already one order of magnitude higher than that before the MRI is triggered. During the TI fully active stage (solid curve), the relatively low surface density region within 0.3 AU corresponds to the TI high stage with T_c > 10⁴ K. The density fluctuations around 0.3 AU result from convection at the transition between the TI high and low states as seen in Figures 1 and 3. We discuss the details of convection in Section 4.3

Figure 5 shows the effective temperature distributions at the same stages as above. When the MRI is triggered, the effective temperature increases dramatically with the increasing accretion rate. Eventually, the whole inner disk almost accretes at a constant rate and the effective temperature is similar to that of a steady accretion solution at $M\dot{M}=1.8\times 10^{-4}\ M_{\odot }^{2}\ {\rm yr}^{-1}$ (smooth solid curve). The discrepancy at small radii may be caused by the outflow inner boundary condition of the two-dimensional simulation, which means all the velocity gradients relative to the coordinates (R, θ) are 0.

3.3. Resolution Study

Figure 6 shows the mass accretion rates during an outburst for models with four different numerical resolutions. In general, the vertical linear resolution must be higher than the radial linear resolution to resolve the vertical structure of the disk, particularly the exponential density drop in the disk photosphere. With N_R = 320 grid cells and N_θ = 112 grid cells, every grid is five times longer radially than vertically.

Zoom In Zoom Out Reset image size

The simulations have a similar maximum mass accretion rate and outburst timescale as long as N_R ≳ 100 grid cells. This suggests that the models are, crudely speaking, converged.

Notice that models with N_R ⩾ 324 grid cells have complex substructure during outburst. This is caused by convection in the TI high-state region. Lower resolution models barely resolve the convective eddies and thus damp away this substructure.

3.4. Boundary Effects

We have also explored the effect of varying the inner boundary radius R_in on the models. Figure 7 shows $\dot{M}(t)$ for a set of models with 0.06 AU ⩽ R_in ⩽ 0.3 AU, with one two-sided (not equatorially symmetric) model at R_in = 0.06 AU. All models have N_R = 320 and N_θ = 112, except the two-sided model, which has N_θ = 224.

Zoom In Zoom Out Reset image size

The outburst amplitude and duration are nearly independent of R_in as long as R_in is far enough in that the TI can be triggered; the TI is not triggered when R_in = 0.3 AU, and so the outburst is slightly stronger and shorter. One implication of the lack of sensitivity of the outburst profile to R_in is that, while the TI can sharpen the outburst slightly, the outburst timescale and maximum mass are largely set in the outer disk, where the MRI becomes active, how much mass is accumulated during the quiescent stage, and so on.

When the TI is triggered, $\dot{M}(t)$ exhibits strong short-timescale variations (panels (b–d) of Figure 7). This is caused by convection during the transition between the TI high and low states (Figure 8); the large temperature gradients at the transition make it highly convectively unstable by the Schwarzschild criterion (see Section 4 for details).

Zoom In Zoom Out Reset image size

Convection at the high state–low state transition radius changes the character of the $\dot{M}(t)$ variability in two-dimensional models compared to one-dimensional (vertically averaged) models. In one-dimensional evolution models flickering at any radius will propagate through the entire inner disk, making the entire inner disk flicker between the low and high states. Thus, the mass accretion rate varies violently during one-dimensional evolutions. To avoid this, Bell & Lin (1994) used an α prescription designed to make α increase sharply from the low to the high states. In two-dimensional models, on the other hand, flickering propagates both radially and vertically. For example, consider a region near the midplane that wants to jump to the high state will soon experience convective cooling and return to the low state. This acts to localize and dampen low state–high state flickering.

For models with R_in = 0.06 AU (panel (c) in Figure 7), $\dot{M}(t)$ varies less violently than in models with R_in = 0.1 AU because the inner edge is far from the convection associated with the low state–high state transition (Figure 8). As the resolution increases, however, the R_in = 0.06 AU model becomes more and more variable. These models, and especially the convection zones, are not converged even at our highest resolution. This is not a major concern, as the fine structure of convection in our axisymmetric, phenomenological α viscosity model will not be the same as the fine structure of convection in a three-dimensional, MHD turbulent disk; further resolving this structure would not improve the model.

One aspect of the details of convection is worth commenting on. Panel (d) in Figure 7 shows a full-disk model, which exhibits less variability than the half-disk models. Convection in the full-disk models can penetrate, rather than reflecting back from, the midplane (Figure 8). This feature of disk convection is consistent with linear theory (Ruden et al. 1988), which shows that the most unstable mode has no nodes in v_z. Because the maximum mass accretion rates and the outburst duration times are similar for half-disk and full-disk models, most of our simulations have been carried with one-sided models to save computation time. However, we need to use the full-disk models to study convection during the outburst.

Although our two-dimensional simulations exhibit complex time-dependent behavior, the vertical structure of the simulated disk can be well fitted by a simple analytical calculation plus some basic assumptions.

Assuming local energy dissipation and with a given effective temperature and viscosity parameter α, the disk structure can be defined by the following equations. Assuming all the viscosity-generated energy radiates vertically, we have

$$\begin{equation} \frac{9}{4}\nu _{V}\rho \Omega ^{2}=\frac{{\it d}}{{\it d} z}\left(F\right) \,, \end{equation} \tag{ 29 }$$

where the vertical flux F is

$$\begin{equation} F=\frac{4\sigma }{3}\frac{{\it d}T_{\tau }^{4}}{{\it d}{\tau }}\,, \end{equation} \tag{ 30 }$$

the viscosity is

$$\begin{equation} \nu _{V}=\alpha c_s^2/\Omega \,, \end{equation} \tag{ 31 }$$

and

$$\begin{equation} {\it d}\tau =-\rho \kappa (T,p){\it d}z \,\,. \end{equation} \tag{ 32 }$$

With hydrostatic equilibrium,

$$\begin{equation} \frac{{\it d}p}{{\it d}z}=-\frac{GM_{*}\rho z}{R^{3}} \,, \end{equation} \tag{ 33 }$$

and the equation of state,

$$\begin{equation} p=\frac{k_{B}}{\mu }\rho T \,, \end{equation} \tag{ 34 }$$

ρ, T, and p are functions of Z. Equation (29) can be integrated toward the midplane, given T_eff = T(τ = 2/3) as the boundary condition. The integration is halted at z_f, where the total viscous heating is equal to σT⁴_eff. In general, z_f ≠ 0, so we change the initial conditions and repeat the integration until z_f = 0, in which case we have the self-consistent vertical structure solution. This procedure is identical to that described in Zhu et al. (2009); we will call these local vertical multizone (LVMZ) models.

Figure 9 shows the vertical structure at 0.1 AU from the two-dimensional simulation and that from the LVMZ with the same central temperature. Although the two-dimensional model cannot resolve the region where T ∼ 10⁴ K, the two-dimensional model otherwise agrees with the LVMZ very well, which means that at 0.1 AU the local treatment is a good approximation even during the TI high state, and that radial energy transport is not important. We find that, in general, the simulation follows the "S" curve calculated by the LVMZ very well during the whole TI activation process (see Appendix B).

Zoom In Zoom Out Reset image size

In the lower left panel of Figure 9, the pressure plateau around τ ∼ 10 is caused by the rapid increase in opacity with increasing temperature near ∼10⁴ K as H ionizes (see the appendix in Zhu et al. 2009). This means that the optical depth can change significantly, while the pressure remains nearly the same.

The large vertical temperature gradient causes this region to become convectively unstable. As shown in the lower right panel, the dense gas is on top of the thin gas at τ < 100, which is clearly unstable against convection or Rayleigh-Taylor instability. Although this convective region is only at the surface and barely resolved at 0.1 AU, it extends to the whole disk at 0.2 AU, where convective patterns even penetrate the midplane (Figure 10).

Zoom In Zoom Out Reset image size

Convection occurs wherever the Schwarzschild criterion is violated: ∇_S ≡ dlog T/dlog P > (γ − 1)/γ = 0.4 in our case (γ = 5/3). For an optically thick disk, if the opacity can be described as κ = κ₀T^βP (Bell & Lin 1994; Zhu et al. 2009), Rafikov (2007) derives ∇_S > (1 + )/(4 − β). Thus, the disk is convectively unstable when its temperature is in the range with the opacity satisfying (1 + )/(4 − β)>(γ − 1)/γ = 0.4.

Before the outburst, the disk is convectively stable. The active layer has a lower temperature than the dust sublimation temperature. Thus, (1 + )/(4 − β) ≃ 0.3, using the Bell & Lin (1994) opacities. The numerical simulation also shows ∇_S is slightly smaller than (γ − 1)/γ = 0.4, which implies it is convectively stable.⁷ The dead zone is isothermal without any energy generation, which means it is also convectively stable with ∇_S∼ 0.

During the TI high state, however, the disk exhibits a variety of convective features at R < 0.6 AU, where T > 2000 K. The inner disk consists of three distinct convective regions.

(1) At 0.35 AU < R < 0.6 AU, the disk is hot enough that the dust has sublimated and the gas is in molecular form (2000 K<T < 5000 K) with TiO and H₂O opacities dominating the opacity. On the basis of the opacity given by Bell & Lin (1994), (1 + )/(4 − β)∼1 so that the disk is convectively unstable. Our numerical simulations have also shown that ∇_S > 0.4 in these regions (changing from 0.4 at R ∼ 0.6 AU to 1 at R ∼ 1 AU). However, convection turns out to be inefficient in these regions, with convective speeds far less than the sound speeds (lower right panel of Figure 10).

(2) At 0.15 AU < R < 0.35 AU, the midplane is hot enough that hydrogen begins to be ionized. The midplane is thermally unstable, and (1 + )/(4 − β) even becomes negative. The disk heats up until the midplane temperature >2 × 10⁴K and it becomes thermally stable. The surface, however, is still far cooler than 2 × 10⁴K, with the effective temperature ∼2000 K. At some height in the disk, there is a transition where the opacity depends steeply on the temperature and convection is strong. At these radii, the transition is close to the midplane and strong convection stirs the entire disk. The largest scale eddies penetrate the midplane.

Figure 11 shows a velocity map of the convective eddies in region (2). We study them by inserting trace particles. These particles rise up with the hot "red" fingers which are in the TI high state. When they approach the surface, they enter the region which is cool, dense, and in the TI low state. Then they fall with the cool "green" fingers to the midplane where they will be heated up again. The convective velocity can be estimated by assuming the particles in the red fingers are accelerated by the buoyancy force from the surrounding dense green fingers. Consider the red fingers are in TI high state with T ∼ 20,000 K and ρ_r, while the green fingers are in TI low state with T ∼ 2000 K and ρ_g ∼ 10ρ_r. Thus, in the rest frame of the green fingers

Zoom In Zoom Out Reset image size

$$\begin{equation} \nu_{r}\sim \left(\frac{H_{e}a}{2}\right)^{1/2}\sim \left(\frac{H_{e} g_{h}}{2}\frac{\rho _{g}-\rho _{r}}{\rho _{r}}\right)^{1/2}\sim\! \left(4.5 H_{e} g_{h}\right)^{1/2}, \end{equation} \tag{ 35 }$$

where H_e is height of the eddy, g_h = GM_*/R³× H/2. While the velocity of the red fingers in the rest, frame of the green fingers is

$$\begin{equation} \nu_{g}\sim \left(\frac{H_{e} g_{h}}{2}\frac{\rho _{g}-\rho _{r}}{\rho _{g}}\right)\sim \left(0.5 H g_{h}\right)\,. \end{equation} \tag{ 36 }$$

The real eddy velocity should between v_g and v_r. If H/R ∼ 0.1, H_e = 0.5 H and R = 0.2 AU, we derive v_g ∼ 2.3 km s⁻¹ and v_r ∼ 7.1 km s⁻¹. Thus, the average eddy velocity should be around 4 km s⁻¹. The maximum velocity should not be larger than twice the average velocity. The lower left panel shows the velocity of these eddies, whose velocity is around 4 km s⁻¹. The maximum velocity takes place at the surface. Although convection is subsonic at the midplane, it is supersonic at the surface (lower right panel in Figure 10).

(3) At R < 0.15 AU, most of the region is already in the TI high state with T > 2 × 10⁴ K. Thus, (1 + γ)/(4 − β) ∼ 0.3, which is similar to the numerical simulation. The surface region where the high state joins to the low state is, however, still convectively unstable (refer to Figure 9 which shows ∇_S is very large at the transition). Unfortunately, our simulation barely resolves this region.

Generally, convection in our simulation is strongly correlated with TI (which is not surprising, since both are related to steep gradients in the opacity). In the outer region where the TI is inactive, convection is weak and inefficient. When the TI is active, the convection patterns depend on the vertical height of the point where hydrogen begins to ionize. If this transition is close to the midplane, convection is strong and penetrates the midplane. Otherwise, only the surface is affected by convection.

Using our adopted parameters, our two-dimensional simulations produce outbursts which are similar to FU Orionis outbursts, with comparable maximum mass accretion rates (∼2 × 10⁻⁴ M_☉ yr⁻¹) and decay timescales (∼100 yr). The outbursts begin with the GI transferring mass to the inner disk until it is hot enough to activate the MRI; at that point the disk accretes at a much higher mass accretion rate because of the much higher efficiency of the MRI in transferring angular momentum.

Our outburst scenario is similar to that of Armitage et al. (2001). However, their outbursts were at a much lower mass accretion rate, ∼10⁻⁵ M_☉ yr⁻¹, and lasted a longer time, ∼10⁴ yr. As discovered by Book & Hartmann (2005) and explained in Zhu et al. (2009), their weaker and slower outbursts are due to their lower MRI trigger temperatures (T_M ∼ 800 K) and smaller α_M = 0.01. Specifically, with a lower T_M, the MRI is triggered at a larger radius with a lower surface density, resulting in less mass accumulation in the inner disk. Combined with a smaller α, this results in a lower mass loss rate and a longer viscous timescale.

Our choice of T_M ∼ 1200 K is close to the temperature range at which dust is expected to sublimate. This suggests that the mechanism for turn-on of the MRI thermally may be controlled by the elimination of the small dust grains that effectively absorb electrons and thus quench the MRI.

Our radiative transfer modeling of FU Ori (Zhu et al. 2007) suggested that to fit the Spitzer Space Telescope IRS spectrum, the rapidly accreting, hot inner disk must extend out to ∼1 AU, inconsistent with a pure TI model (Bell & Lin 1994). In our two-dimensional results, we find MRI triggering at about 2 AU, much more consistent with the Spitzer data and our previous analysis in Zhu et al. (2009). As FU Ori has a smaller estimated central star mass (0.3 M_☉) than the 1 M_☉ used in our simulations, its hot inner disk should be slightly smaller than derived here, but still in agreement with observations.

Non-Keplerian rotation has not been observed from optical to near-IR (5μm) lines (Zhu et al. 2009) since these lines come from inner disk within 1 AU where the disk is not massive. However, due to the fact that the disk is massive and gravitationally unstable at large radii (>2 AU), the disk has slightly sub-Keplerian rotation there and this effect may be observable in the far-IR or submm (Lodato & Bertin 2003).

We cannot calculate meaningful optical (B magnitude) light curves for comparison with the historical observations of FU Ori objects because our inner disk radius is too large; this results in our maximum effective temperatures being too low in comparison with observations. (This constraint was due to the very long computation times required during outburst; moving the inner radius inward greatly lengthens the computing time.) We therefore computed a 1 μm light curve as a prediction to compare with future observations (Figure 12), assuming that the disk radiates like a blackbody at each radius's effective temperature; the lower panel shows the disk's surface and central temperature at 0.06 AU. The outburst has rise time ∼10 yr, which is similar to V 1515 Cyg but longer than FU Ori. The TI is triggered on a much shorter timescale, as the lower panel shows the disk's central temperature rises sharply within one year. However, since the mass at the inner disk is still low, the mass accretion rate and the effective temperature cannot increase significantly until the MRI transfers more mass from the outer disk to the inner disk on a viscous timescale. Overall, our simulation agrees with V 1515 Cyg reasonably well but has a longer rise time than FU Ori and V 1057 Cyg. This might be improved by using a smaller inner radius, a consistent central mass, and a self-consistent boundary condition (Kley & Lin 1999).

Zoom In Zoom Out Reset image size

Our simulation suggests that the short-timescale variations appearing in the light curves of FU Ori (Kenyon et al. 2000) might be caused by convection in the high state–low state transition region. The observations show ΔM ≃ 0.035 and timescales of less than 1 day. The amplitude of the variability in the two-dimensional models is ΔM ∼ 0.2 (see Figure 12). In three dimensions, the amplitude of the variability will be reduced because there will be several convective eddies at each radius and averaging over these will reduce the variability. If we assume that the azimuthal extent of a convective eddy is ∼H, then there will be m = 2πR/H convective eddies and the variability will be reduced by a factor of m^1/2 ∼ 5 if H/R ∼ 0.2, suggesting a variability of ΔM ∼ 0.04. This agreement of the data is encouraging, but it is likely that the details of convection in a magnetized, three-dimensional disk will be different from that found in our viscous, axisymmetric models. The modest success of our model, however, might motivate the development of energetically self-consistent three-dimensional MHD models of the inner disk of FU Ori.

The timescale of variability in our two-dimensional models (months) is much longer than the observed 1 day. The timescale of variability might be shortened by (1) increasing resolution, since the timescale of variability in our models increases with resolution; (2) switching off all components of stress in Equation (4) except the r–Φ component, since the other components act to retard and smooth convection; (3) replacing the viscosity model used here by MHD turbulence, since the viscosity tends to suppress small scale structure; (4) three-dimensional effects, since one might expect to see variability on timescales of 2π/(mΩ) ∼ 1 day, using R = 0.2 AU and our estimate for m from above. Overall, although convective behavior revealed by our current simulations are very preliminary (they depend on the simulation resolution, opacity, details of the turbulence, etc.), it seems reasonable that convection accounts for the short-timescale variability of FU Orionis objects.

Hartmann et al. (2004) found that some mildly supersonic turbulence (∼2 c_s) was needed to explain the width of ¹²CO first-overtone lines of FU Ori. High-resolution spectral modeling suggests that CO lines broadened with ∼4 km s⁻¹ turbulent velocity fit the observations well (Zhu et al. 2007). Hartmann et al. (2004) proposed that the MRI could be a mechanism to drive the supersonic turbulence, while here we suggest that convection could also drive the turbulence. Convection is subsonic at the midplane in the TI high state, but as the convective eddies travel to the surface the temperature drops and convection becomes supersonic (see the lower right panel of Figure 10). Numerical simulations (lower left panel of Figure 10) and the simple analytical estimate (Equation (36)) suggest the 5 km s⁻¹ is a typical convective velocity.

We also studied the disk central temperature at the stage when the MRI is activated and then moves inward. The transition between the TI active and inactive region is sharp and occurs over a radial interval comparable to the disk-scale height. The fronts travel at 2.7 km s⁻¹, which is almost the sound speed at T = 1000 K. This thermal front velocity is much higher than αc_s predicted by Lin et al. (1985), and we will discuss this phenomenon in detail in our next parameter study paper. Notice that in our model the stress is assumed to respond instantaneously to a change in the temperature caused by the radial energy diffusion; this may indeed be the case if magnetic field mixes into the newly heated gas from the hot side of the transition front. If the magnetic field takes longer than ∼Ω⁻¹ to build up, then the front will travel more slowly.

Radiation pressure starts to become significant during the TI high state, with P_r/P_t∼ 0.3 in the innermost regions. We have conducted some experiments including radiation pressure using the diffusion approximation and find generally similar behavior to what has been shown in this paper, though the disk is slightly thicker with the inclusion of radiation pressure. Again, we note that our inner radius is considerably larger than those of T Tauri stars, so radiation pressure might be more important as accreting material joins the central star.

We have studied the time evolution of FU Orionis outbursts in a layered disk scenario, using an axisymmetric, two-dimensional viscous hydrodynamics code with radiative transfer. Both the MRI and GI are considered as angular momentum transport mechanisms and are modeled by a phenomenological α prescription. Besides the surface of the disk ionized by cosmic rays and X-rays, the MRI is also considered to be active as long as the disk's temperature is above a critical temperature T_M so that thermal ionization can provide enough ions to couple the plasma to the magnetic field.

Starting with a massive disk equivalent to a constant α disk with $\dot{M}=10^{-4}$ M_☉ yr⁻¹ and α = 0.1, several outbursts appear during our simulation. These outbursts begin with the GI transferring mass to the inner disk until the inner disk gradually becomes gravitationally unstable. Eventually, the disk at 2 AU is hot enough to trigger the MRI and then the MRI active region expands throughout the whole inner disk. With the increasing temperature since the MRI is transferring more mass to the inner disk, the TI is triggered and a "high state" expands outward to 0.3 AU. The outburst state lasts for hundreds of years and drains the inner disk; the disk then returns to the low state.

Convection appears at the inner disk during the TI high state. In the outer region (0.35 < R < 0.6 AU), where the TI is inactive, convection is weak and inefficient. At 0.15 < R < 0.35 AU during outburst, there is a hydrogen ionization zone at z ∼ H, and strong convection is present. The convection even penetrates the midplane to form antisymmetric convective rolls. In the innermost region, most of the disk is hot, the hydrogen ionization zone is close to the photosphere, and weak convection is confined to the disk surface.

Our model shows maximum accretion rates and outburst duration timescales similar to those of the known FU Orionis objects. Unlike the pure TI theory, our model produces an AU-scale inner disk with high mass accretion rate, in agreement with the observational constraints.

With some simple assumptions our simulations have produced FU Orionis-type outbursts. They also pose some interesting questions for future study. In a later paper, we will compare our two-dimensional results to one-dimensional simulations; the latter will allow us to consider longer disk evolutionary timescales.

A dynamical GI simulation at AU scale, with both self-gravity and radiation, will be important to understand angular momentum transport by GI and the thermal structure of the inner disk. However, considering the large Keplerian velocity of the inner disk, a self-consistent radiation hydrodynamic scheme will be needed. Our simulations have also shown that convection is important in FU Orionis objects due to the large vertical temperature gradient at the hydrogen ionization zone. A fully three-dimensional simulation with both MHD and radiation, together with a self-consistent treatment of the equation of state, will be essential to understanding the role of convection. Also the very inner disk where the MRI is thermally activated by the irradiation and the boundary layer may also play a role in the disk evolution.

Z.Z. wants to thank Dr. R. Durisen and K. Cai for the inspiring discussion about their radiative transport scheme. This work was supported in part by NASA grant NNX08A139G, by the University of Michigan, and by a Sony Faculty Fellowship, a Richard and Margaret Romano Professorial Scholarship, and a University Scholar appointment to CG. JCM was supported by NASA's Chandra Fellowship PF7-80048.

To test the radiative transfer scheme (and particularly the boundary condition at the disk photosphere), we have constructed a disk model in which the viscous heating is turned off when τ < 100. This makes the vertical flux constant for τ < 100 and allows us to test our scheme against a gray atmosphere model. The disk structure shown here is that calculated with Method 1 as discussed in Section 2.2. The structure calculated with Method 2 is even in a better agreement than Method 1, since it directly uses the gray atmosphere structure as the boundary condition.

The left panel in Figure 13 shows the flux through the atmosphere at different heights for one radius, while the dotted line is F = σT⁴_eff with T_eff derived from the simulation. The right panel shows the disk temperature structure at this radius from the two-dimensional simulation (solid line) and the theoretical gray atmosphere temperature structure (T⁴ = 3/4 T⁴_eff(τ+2/3)) (dotted line). The vertical structure in the two-dimensional simulation agrees well with the gray atmosphere structure and improves on the original scheme of Cai et al. (2008) by eliminating a temperature jump at the photosphere.

Zoom In Zoom Out Reset image size

The TI model was developed to account for outbursts in dwarf novae systems (Faulkner et al. 1983). But it also has features that make it attractive for explaining FU Ori outbursts (Bell & Lin 1994). Although in this paper we find that MRI triggering is mainly responsible for initiating FU Orionis outbursts, the TI still exists and can be triggered during the outburst. With the increasing mass accretion rate by the MRI operating, the innermost disk is heated up to hydrogen ionization and the disk becomes thermally unstable. Finally, TI sets in and the innermost disk enter a high (fully ionized) state.

The basic idea behind the TI can be illustrated by the disk Σ − T_c relationship ("S" curve, Figure14). The disk is in thermal equilibrium on this curve. To the left of the curve cooling exceeds heating and to the right heating exceeds cooling. As the disk central temperature reaches 3000 K the disk transitions from the low-state lower branch to the high-state (T_c ∼ 2 × 10⁴ K) upper branch; it cannot sit on the unstable middle branch.

Zoom In Zoom Out Reset image size

B.1. "S" Curve of Two-dimensional Models

B.2. "S" Curve of One-dimensional Models

In one zone model, where the midplane opacity is used at all z and the total surface density replaces ρ(z), the "S" curve can be obtained via

$$\begin{equation} F_{{\rm vis}}=F_{{\rm rad}} \,. \end{equation} \tag{ 37 }$$

Thus,

$$\begin{equation} \frac{9}{4}\alpha \Omega \Sigma c_{s}^{2}=2\sigma T_{{\rm eff}}^{4}\,, \end{equation} \tag{ 38 }$$

where Σ is the surface density and c_s is the sound speed at the local central disk. For an optically thick disk, we have T⁴_c = $\frac{3}{8}\tau _{R}$T⁴_eff, where τ_R = κ_RΣ and κ_R is the Rosseland mean opacity (Hubeny 1990). Then,

$$\begin{equation} T_{c}^{3}=\frac{27}{64\mu \sigma }\alpha \Omega \Sigma ^{2} {\cal {R}}_{c} \kappa _{R} \,, \end{equation} \tag{ 39 }$$

where $\cal {R}_{\rm c}$ is the gas constant and μ is the mean molecular weight. This T_c–Σ relationship is shown in Figure 14, which differs significantly from the LVMZ "S" curve.

These "S" curves again make the point that the observed large outer radius, 0.5–1 AU, of the high accretion rate portion of the FU Ori disk poses difficulties for the pure TI model. Even in the situation explored by Lodato & Clarke (2004), in which the TI is triggered by a massive planet in the inner disk, the outer radius of the high state is still the same as in Bell & Lin (1994), ∼ 20 R_☉. The relatively high temperatures (2 − 4 × 10³ K) required to trigger the TI (due to the ionization of hydrogen) are difficult to achieve at large radii; they require large surface densities which in turn produce large optical depths, trapping the radiation internally and making the central temperature much higher than the surface (effective) temperature.