CONVECTION CAUSES ENHANCED MAGNETIC TURBULENCE IN ACCRETION DISKS IN OUTBURST

The basic equations for our radiation MHD simulations are

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol {v}) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial (\rho \boldsymbol {v})}{\partial t} + \nabla \cdot (\rho \boldsymbol {v}\boldsymbol {v}) = -\nabla p + \frac{1}{4\pi }(\nabla \times \boldsymbol {B})\times \boldsymbol {B} + \frac{{\kappa }_{\rm R}\rho }{c}\boldsymbol {F}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot (e\boldsymbol {v}) = -(\nabla \cdot \boldsymbol {v})p - \left(4\pi B(T) - cE\right){\kappa }_{\rm P}\rho, \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial E}{\partial t} + \nabla \cdot (E\boldsymbol {v}) = - \nabla \boldsymbol {v}:\mathsf {P} + \left(4\pi B(T) - cE\right){\kappa }_{\rm P}\rho - \nabla \cdot \boldsymbol {F}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{\partial \boldsymbol {B}}{\partial t} - \nabla \times \left(\boldsymbol {v}\times \boldsymbol {B}\right) = 0, \end{equation} \tag{ 5 }$$

where ρ is the gas density, e the gas internal energy, p the gas pressure, T the gas temperature, E the radiation energy density, $\mathsf {P}$ the radiation pressure tensor, $\boldsymbol {F}$ the radiation energy flux, $\boldsymbol {v}$ the velocity field vector, $\boldsymbol {B}$ the magnetic field vector (in CGS emu units), B(T) = σ_BT⁴/π the Planck function (σ_B, the Stefan–Boltzmann constant), and c the speed of light. We use a flux-limited diffusion approximation of the radiative transfer, where the energy flux $\boldsymbol {F}$ and pressure tensor $\mathsf {P}$ are related to the energy density E as $\boldsymbol {F} = -(c\lambda (R)/\kappa _{\rm R}\rho)\nabla E$ and $\mathsf {P} = \mathsf {f}(R)E$. Here λ(R) ≡ (2 + R)/(6 + 3R + R²) is a flux limiter with R ≡ |∇E|/(κ_RρE), and $\mathsf {f}(R) \equiv (1/2)(1-f(R))\mathsf {I} + (1/2)(3-f(R))\boldsymbol {n}\boldsymbol {n}$ is the Eddington tensor with f(R) ≡ λ(R) + λ(R)²R² and $\boldsymbol {n}\equiv \nabla E/|\nabla E|$ (Turner & Stone 2001).

The EOSs, p = p(ρ, e/ρ) and T = T(ρ, e/ρ), are computed from the Saha equations assuming ionization equilibrium with solar abundances (Figure 1). Our method is equivalent to that described in Appendix A in Tomida et al. (2013) with the following exceptions: (1) we consider metals and H$_2^+$ in addition to the species employed in their method and (2) we allow the orthohydrogen to parahydrogen ratio to be determined by thermal equilibrium rather than assuming the high temperature fixed ratio of 3:1. (Neither of these changes make significant differences.)

Zoom In Zoom Out Reset image size

The Planck-mean opacity κ_P and Rosseland-mean opacity κ_R are given as a function of density ρ and gas temperature T (Figure 2). We have combined three published opacity tables, Semenov et al. (2003), Ferguson et al. (2005), and the Opacity Project (OP; Seaton 2005), to cover the relevant temperature range for our purpose (10³ < T < 10⁶ K). Different opacity tables are connected (with linear interpolation) where they seem mostly consistent; specifically, the Semenov et al. and Ferguson et al. opacities are connected at the dust sublimation temperatures while the Ferguson et al. and OP opacities are connected at T = 10^3.7 K. The combined opacity tables have upper and lower bounds, beyond which the opacities are extended using a zero-gradient extrapolation, which is not a bad approximation for T > 10³ K (see Figure 1 in Malygin et al. 2013).

Zoom In Zoom Out Reset image size

In the shearing box approximation, a local patch of an accretion disk is modeled as a co-rotating Cartesian frame (x, y, z) with linearized Keplerian shear flow −(3/2)Ωx, where the x, y, and z directions correspond to the radial, azimuthal, and vertical directions, respectively (Hawley et al. 1995). In this approximation, the inertial force terms $-2\Omega \hat{\boldsymbol {z}}\times \boldsymbol {v} +3\Omega ^2x\hat{\boldsymbol {x}} - \Omega ^2z\hat{\boldsymbol {z}}$ (Coriolis force + tidal force + vertical component of gravitational force, respectively) are added to the equation of motion (2), where $\hat{\boldsymbol {x}}$ and $\hat{\boldsymbol {z}}$ are the unit vectors in the x and z direction, respectively. Shearing-periodic, periodic, and outflow boundary conditions are used for the x, y, and z boundaries of the box, respectively (Hirose et al. 2006).

The radiation MHD equations are solved time-explicitly by ZEUS using the Method of Characteristics–Constrained Transport algorithm except for the radiation–gas energy exchange terms ±(4πB − cE)κ_Pρ and the radiative diffusion term $-\nabla \cdot \boldsymbol {F}$ whose time scales are much shorter than the MHD time scale (Stone & Norman 1992a, 1992b; Turner & Stone 2001). Those terms are coupled and solved time-implicitly using Newton–Raphson iteration and the multi-grid method with a Gauss–Seidel smoother (Tomida et al. 2013). The gas temperature T used in evaluating the mean opacities is fixed to that in the previous time step to linearize the radiative diffusion equation.

We assume no explicit resistivity and viscosity in the basic equations, and thus turbulent dissipation occurs through the sub-grid numerical dissipation. The numerically dissipated energy is captured in the form of additional internal energy in the gas, effectively resulting in an additional term Q_diss in the gas energy equation (3). To accomplish this, the original ZEUS algorithm is modified so as to conserve total energy (Turner et al. 2003; Hirose et al. 2006).

During the simulations, we employ a density floor of 10⁻⁶ of the initial midplane density to avoid very small time steps. We also employ a small internal energy floor for numerical stability (Hirose et al. 2006). The total artificial energy injection rate associated with these floors and numerical errors in the implicit solver is generally less than ∼1% of the cooling/heating rates of the final steady state.

The initial conditions within the disk photosphere are determined as follows: first, the vertical profiles of mass density ρ(z), pressure p(z) and temperature T(z), as well as the initial surface density Σ₀ and the photosphere height h₀, are determined from an α model described in the Appendix, by choosing an alpha α₀ and an effective temperature T_eff0. Then, the internal energy density is calculated via the EOS e(z) = e(p(z), T(z)), and the radiation energy density is calculated by E(z) = 4σ_BT⁴(z)/c, assuming that the gas and radiation are in thermal equilibrium with each other. Above the disk photosphere z = h₀, we assume a uniform, low-density atmosphere with ρ(z > h₀) = ρ(0) × 10⁻⁶, e(z > h₀) = e(h₀) × 10⁻⁶, and E(z > h₀) = E(h₀) × 10⁻⁶.

The initial velocity field is the linearized Keplerian shear flow $\boldsymbol {v} = (0,-3/2\Omega x,0)$, whose x and z components are perturbed by random noise of 0.5% of the local sound velocity. The initial configuration of the magnetic field is a twisted azimuthal flux tube (with net azimuthal flux, but zero net vertical flux), which is placed at the center of the simulation box and is confined within the photospheres. (Note that the volume-integrated vertical magnetic flux is conserved to be zero during the simulation while the volume-integrated azimuthal magnetic flux is not conserved since it can escape from the box through the vertical boundaries.) The field strength in the tube is uniform and the ratio of the poloidal field to the total field is 0.25 at maximum. The field strength is specified by the parameter β₀, the ratio of thermal pressure to magnetic pressure at the center of the flux tube. The results do not depend on β₀ or the shape of the flux tube so long as the initial development of the MRI is well resolved.

Parameters of the simulations are listed in Table 1. As described above, the initial conditions are specified by two parameters, surface density Σ₀ (or equivalently α₀) and effective temperature T_eff0. The table also lists time-averaged surface density $\bar{\Sigma }$, effective temperature $\bar{T}_{\rm eff}$, and α in the final steady state (see Section 3 for the definitions of these quantities).

Table 1. List of Runs

Run	Σ0	Teff0	α0	h0	β0	$\bar{\Sigma }$	$\bar{T}_{\rm eff}$	α	τtot	Nx	Ny	Nz	Lx/h0	Ly/h0	Lz/h0	Lz/hp	tth	t1	t2
Upper Branch Solutions
ws0430F	2983	14454	0.0097	1.81E+09	100	2900	20617	0.0401	37931	32	64	256	0.500	2.000	4.000	8.37	6.53	50	150
ws0429F	1075	10000	0.0098	1.41E+09	100	1030	13352	0.0332	24232	32	64	256	0.500	2.000	4.000	8.54	9.49	90	190
wa0429H	1075	10000	0.0098	1.41E+09	100	1070	12570	0.0279	29083	48	96	384	0.600	2.400	4.800	10.4	11.5	90	190
wv0429W	1075	10000	0.0098	1.41E+09	100	1052	13067	0.0311	26674	64	128	256	1.000	4.000	4.000	8.58	10.2	20	120
ws0439F	555	7943	0.0099	1.12E+09	100	540	10649	0.0325	20337	32	64	256	0.500	2.000	4.000	7.79	10.9	50	150
ws0469F	402	7079	0.0098	9.74E+08	100	391	9514	0.0323	19742	32	64	256	0.560	2.240	4.480	8.14	11.6	50	150
wa0469H	402	7079	0.0098	9.74E+08	100	400	9463	0.0318	21031	48	96	384	0.672	2.688	5.376	9.75	11.7	50	150
wv0469W	402	7079	0.0098	9.74E+08	100	385	9386	0.0328	22161	64	128	256	1.120	4.480	4.480	8.23	11.5	50	150
ws0468C	397	5011	0.0031	7.40E+08	100	386	9594	0.0375	22518	32	64	256	0.750	3.000	6.000	8.37	10.2	50	150
ws0491F	342	12882	0.0988	1.27E+09	100	332	8992	0.0347	22712	32	64	256	0.500	2.000	4.000	9.77	11.3	50	150
ws0470F	281	12022	0.0991	1.21E+09	100	273	9059	0.0452	17420	32	64	256	0.500	2.000	4.000	9.81	8.74	50	150
ws0472C	248	4073	0.0031	5.15E+08	100	239	8492	0.0459	23321	32	64	256	1.140	4.560	9.120	10.4	9.06	50	150
ws0471F	247	11481	0.0988	1.18E+09	100	238	8853	0.0522	18117	32	64	256	0.500	2.000	4.000	10.0	7.84	50	150
wa0471H	247	11481	0.0988	1.18E+09	100	243	8872	0.0496	29032	48	96	384	0.600	2.400	4.800	12.0	8.06	20	120
wv0471W	247	11481	0.0988	1.18E+09	100	235	8766	0.0546	23992	64	128	256	1.000	4.000	4.000	10.3	7.69	50	150
ws0492F	217	10964	0.0986	1.14E+09	100	211	8366	0.0548	23980	32	64	256	0.500	2.000	4.000	10.5	7.90	50	150
ws0425F	204	10715	0.0985	1.11E+09	100	197	8582	0.0668	20029	32	64	256	0.500	2.000	4.000	10.1	6.48	50	150
ws0427F	191	10471	0.0981	1.09E+09	100	185	8482	0.0646	16517	32	64	256	0.500	2.000	4.000	10.0	6.66	50	150
ws0437F	179	10232	0.0975	1.07E+09	100	174	8283	0.0651	16791	32	64	256	0.500	2.000	4.000	10.1	6.64	50	150
wa0437H	179	10232	0.0975	1.07E+09	100	177	8908	0.0821	14401	48	96	384	0.600	2.400	4.800	11.5	5.17	50	150
wv0437W	179	10232	0.0975	1.07E+09	100	168	8333	0.0795	23234	64	128	256	1.000	4.000	4.000	10.6	5.79	50	150
ws0433F	168	10000	0.0978	1.05E+09	100	162	8263	0.0749	15595	32	64	256	0.500	2.000	4.000	10.2	5.99	50	150
ws0436F	158	9772	0.0978	1.03E+09	100	154	8138	0.0747	17608	32	64	256	0.500	2.000	4.000	10.3	6.05	30	130
wt0487F	149	9549	0.0983	1.00E+09	100	143	7914	0.0862	20503	32	64	256	0.500	2.000	4.000	10.6	5.58	50	150
ws0441F	140	9332	0.0981	9.77E+08	100	134	7966	0.0927	17556	32	64	256	0.500	2.000	4.000	10.2	5.12	50	150
ws0494C	140	6606	0.0311	7.07E+08	10	138	7868	0.0919	23865	32	64	256	0.700	2.800	5.600	10.8	5.29	10	90
ws0446F	132	9120	0.0987	9.51E+08	100	127	7490	0.1062	34556	32	64	256	0.500	2.000	4.000	12.0	5.06	50	150
wa0446H	132	9120	0.0987	9.51E+08	100	128	7876	0.1087	23795	48	96	384	0.600	2.400	4.800	13.0	4.66	50	150
wz0446W	132	9120	0.0987	9.51E+08	10	123	7758	0.1210	28978	64	128	256	1.000	4.000	4.000	11.4	4.35	25	125
wt0442F	118	8709	0.0982	8.97E+08	10	113	7522	0.1195	22037	32	64	256	0.500	2.000	4.000	11.4	4.53	50	150
ws0488R	99	8128	0.0982	8.12E+08	10	—	—	—	—	32	64	256	0.500	2.000	4.000	—	—	—	—
Lower branch solutions
ws0467R	373	2511	0.0031	1.24E+08	100	—	—	—	—	32	64	256	1.000	4.000	8.000	—	—	—	—
ws0466F	353	1584	0.0031	8.18E+07	100	348	2675	0.0204	13.2	32	64	256	0.750	3.000	6.000	8.18	17.3	20	100
ws0438F	285	1096	0.0010	1.13E+08	100	275	2714	0.0287	8.8	32	64	256	0.500	2.000	4.000	7.70	12.3	50	150
wt0435F	211	1000	0.0010	1.06E+08	100	191	2271	0.0271	2.6	32	64	256	0.450	1.800	3.600	7.44	10.9	120	220
ws0462F	176	1778	0.0099	6.51E+07	100	174	2283	0.0282	2.4	32	64	256	1.000	4.000	8.000	10.1	9.78	50	150
ws0464C	121	3388	0.0313	1.35E+08	100	117	1969	0.0254	2.7	32	64	256	0.500	2.000	4.000	10.9	10.6	50	150
ws0465F	120	2344	0.0314	5.18E+07	100	116	2051	0.0312	2.3	32	64	256	1.000	4.000	8.000	8.29	8.89	50	150
ws0434F	100	1513	0.0099	5.85E+07	100	93	1885	0.0295	3.2	32	64	256	0.750	3.000	6.000	7.12	9.75	50	150
ws0476F	49	1698	0.0315	5.28E+07	100	45	1560	0.0314	5.4	32	64	256	0.800	3.200	6.400	7.34	9.02	50	150
ws0445F	32	1096	0.0098	8.63E+07	100	31	1515	0.0418	4.4	32	64	256	0.500	2.000	4.000	7.52	6.66	50	150

Notes. The last letters in the names of runs denote the following: F = fiducial run, H = high resolution run, W = wide box run, R = runaway heating/cooling run, and C = run to check dependence on the initial condition. The units of surface densities (Σ₀ and $\bar{\Sigma }$), effective temperatures (T_eff0 and $\bar{T}_{\rm eff}$), height (h₀), and thermal time (t_th) are, respectively, g cm⁻², K, cm, and orbit. L_x, L_y, and L_z are the lengths, and N_x, N_y, and N_z are the numbers of cells, in the x, y, and in z directions, respectively. The time-averaging in diagnostics is done for t₁ < t < t₂ orbits. The pressure scale height of the steady state is computed as h_p ≡ ∫[ < p_thermal > ]dz/2max ([ < p_thermal > ]).

Download table as:  ASCIITypeset image

There are also two numerical parameters: the box size and the number of cells. Because we consider dependence on the surface density in this work, it is desirable that the surface density does not change much during the simulation; therefore we adjust the box height so that the mass loss rate through the top/bottom boundaries is less than about 10% of the initial mass per hundred orbits, while the MRI is kept resolved reasonably near the disk midplane. Such a box is typically ∼10 scale heights tall when measured in terms of the final steady-state pressure scale height h_p (see Table 1). In the fiducial runs, which are the ones discussed in Section 3, the numbers of cells are (32, 64, 256) and the aspect ratio of the box is 1:4:8 in the x, y, and z directions, respectively. See the table for the absolute size of the box. To check numerical convergence, we have also run a wide box version (the box size in both the x and y directions is doubled) and a high resolution version (the resolution is 1.5 times higher and the box length in each direction is 1.2 times larger) for some selected fiducial runs, which we discuss in Section 4.2.

One way to characterize the results is in terms of the effective temperature at the disk photosphere. Figure 3 shows this quantity as a function of surface density for equilibrium states that last over 100 orbits.⁷ The time-averaged effective temperature and surface density are computed as

$$\begin{eqnarray} & {\bar{T}_{\rm eff} \equiv \left[\left(\frac{1}{2\sigma _{\rm B}}\int \langle Q^- \rangle dz\right)^{\frac{1}{4}}\right]}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} & \bar{\Sigma } \equiv \left[\int \left<\rho \right>dz\right],\qquad\ \ \ \ \end{eqnarray} \tag{ 8 }$$

where 〈Q⁻〉 is the total cooling rate,

$$\begin{eqnarray} & \langle Q^- \rangle \equiv \frac{d}{dz}\left<F_z\right> + \frac{d}{dz}\left<(e + E)v_z\right>. \end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

The brackets [ ] here denote time averaging over a selected period of 100 orbits in which the disk is in quasi-steady state and the MRI in the disk is fairly resolved (the period of time averaging in each run is listed in Table 1), and the space integration is done for the full extent of the box in z. As mentioned in the previous section, the surface density can vary due to the mass loss through the vertical boundaries, and thus the time averaged surface density $\bar{\Sigma }$ is typically smaller than the initial surface density Σ₀ by only a few percent, rising to at most 9% (see Table 1).

The DIMs based on the α prescription of the stress generally produce S-shaped thermal equilibrium curves in this plane, as illustrated by the gray curves. Our simulations show that "S-curves" can also arise as a result of MRI turbulence at temperatures near the hydrogen ionization transition as reported by Latter & Papaloizou (2012). There are two major solution branches: the upper hot branch ($\bar{T}_{\rm eff} \gtrsim 8000$ K, $\bar{\Sigma } \gtrsim 100$ g cm⁻²) and the lower cool branch ($\bar{T}_{\rm eff} \lesssim 3000$ K, $\bar{\Sigma } \lesssim 300$ g cm⁻²). For a limited range of surface density ($100 \lesssim \bar{\Sigma } \lesssim 300$ g cm⁻²), there exist two different stable states for a single value of surface density, showing bistability.

The disk is almost fully ionized and optically thick (with total optical depth τ_tot > 10⁴) on the upper branch, but it is almost wholly neutral and much less optically thick (2 < τ_tot < 14) on the lower branch. (See Table 1 for the value of τ_tot in each run.) These features of the upper and lower branches are consistent with the DIMs.

The detailed structures of gas temperature as well as ionization fraction are given in Figure 4 for three selected runs: ws0429F, ws0446F, and ws0465F. (Radiation temperature is almost identical to gas temperature in these cases and thus is not drawn.) The temperature is above 10⁴ K and the ionization fraction is almost unity on the extreme right end of the upper branch (Figure 4(A)). However, as the surface density decreases toward the left edge of the upper branch, the temperature and ionization fraction can fall from values in this range near the midplane down to ∼6000 K and ∼10⁻², respectively, in the atmosphere (Figure 4(B)). We expect that this low ionization in the atmosphere will not significantly affect MRI turbulence near the midplane, where the gas is fully ionized.⁸ On the lower branch, the temperature is typically below 2000 K and the ionization fraction drops to <10⁻⁷ (Figure 4(C)). Here, the vertical temperature gradient is very shallow due to the low optical depth.

Zoom In Zoom Out Reset image size

A new finding here is that α is not constant, as shown by the colors in Figure 3. To be precise, we define α as

$$\begin{eqnarray} &&\quad\ \ \alpha \equiv \frac{[W_{xy}]}{[P_{\rm thermal}]}, \end{eqnarray} \tag{ 10 }$$

where vertically integrated total stress W_xy and thermal pressure P_thermal are defined as

$$\begin{eqnarray} &{W}_{xy}(t) \equiv \int \langle w_{xy} \rangle dz, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \, &{P}_{\rm thermal}(t) \equiv \int \left<p_{\rm thermal}\right>dz. \end{eqnarray} \tag{ 12 }$$

Here p_thermal ≡ p + E/3 and w_xy ≡ −B_xB_y + ρv_xδv_y, where δv_y ≡ v_y + (3/2)Ωx. (We include radiation pressure E/3 in the thermal pressure, but its fraction is at most 0.057 in the largest surface density case on the upper branch.) Changing the definition of α to a time-average of the instantaneous ratio, [W_xy/P_thermal], changes the values by typically a few percent, and always less than 10% for those simulations achieving a quasi-steady state.

Throughout the lower branch and at the right end of the upper branch, α is approximately 0.03, typical of values found previously in local numerical simulations of MRI turbulence that lack net vertical magnetic flux. However, near the left edge of the upper branch (7000 ≲ T_eff ≲ 10, 000 K), α rises to as much as 0.12, increasing as the surface density decreases. This behavior is in contrast with the DIMs, where constant values of α are assumed on each branch: ∼0.1 on the upper branch and ∼0.01 on the lower branch (see Section 4.3). As shown below, vertical profiles of the energy transport reveal that the notable change in α in our simulations is associated with thermal convection.

In Figure 5, time-averaged profiles of radiative heat flux $\bar{F}^-_{\rm rad}(z)$, advective heat flux $\bar{F}^-_{\rm adv}(z)$, and cumulative heating rate $\bar{F}^+_{\rm heat}(z)$ are shown for the cases treated in Figure 4. Here, the heat fluxes are defined as

$$\begin{eqnarray} &&\hspace{-164pt} \bar{F}^-_{\rm rad}(z) \equiv \left[\left<F_z\right>\right], \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\hspace{-138pt} \bar{F}^-_{\rm adv}(z) \equiv \left[\left<(e + E)v_z\right>\right], \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} & \bar{F}^+_{\rm heat}(z)\equiv \int _0^z\left[\left<Q_{\rm diss}\right> + \left<-\left(\nabla \cdot \boldsymbol {v}\right)p\right> + \left<-\nabla \boldsymbol {v}:\mathsf {P}\right>\right]dz. \quad\quad \end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

From the energy equations (3) and (4), the thermal energy balance in a steady state is written as

$$\begin{eqnarray} &&Q_{\rm diss} - \left(\nabla \cdot \boldsymbol {v}\right)p - \nabla \boldsymbol {v}:\mathsf {P} = \nabla \cdot \boldsymbol {F} +\nabla \cdot \left((e + E)\boldsymbol {v}\right), \quad \end{eqnarray} \tag{ 16 }$$

where the left hand side is the heating rate (turbulent dissipation and compressional heating) while the right hand side is the cooling rate (radiative diffusion and advection). Therefore, it is expected in Figure 5 that

$$\begin{eqnarray} &&\bar{F}^+_{\rm heat}(z) = \bar{F}^-_{\rm rad}(z) + \bar{F}^-_{\rm adv}(z), \end{eqnarray} \tag{ 17 }$$

which is a vertically integrated form of Equation (16). Actually, Equation (17) roughly holds in panel (B), and almost exactly holds in panels (A) and (C). (Note that the gray curve ($\bar{F}^+_{\rm heat}$) almost matches the red curve ($\bar{F}^-_{\rm rad}$) in panels (A) and (C))

When radiative diffusion carries the dissipated turbulent energy to the disk surface (Figures 5(A) and (C)), the value of α is typical of MRI turbulence. On the other hand, when advection plays a major role in the vertical heat transport near the disk midplane (Figure 5(B)), α is large. As shown in the next subsection, this advective heat transport is associated with hydrodynamic thermal convection, and it is triggered by high opacity that suppresses heat transport by radiative diffusion. We therefore refer to this advective heat transport as convection from now on. The energy transported upward by convection is deposited at higher altitude, but below the photosphere, and is then carried by radiative diffusion toward the disk surface.

As shown in Figures 6(A) and (B), the left edge of the upper branch, where α is large, is exactly the place where convection carries much of the heat flux. Here, to measure the convection strength, we define the advective fraction as the ratio of the advective energy flux to the total energy flux,

$$\begin{eqnarray} &&f_{\rm adv} \equiv \frac{\int [\langle (e+E)v_z \rangle] {\rm sgn}(z) [\langle p_{\rm thermal} \rangle]dz}{\int \left[\left<\left(e+E\right)v_z\right> + \left<F_z\right>\right]{\rm sgn}(z)\left[\left<p_{\rm thermal}\right>\right]dz}, \end{eqnarray} \tag{ 18 }$$

where the time-averaged pressure [ < p_thermal > ] is used as a weight function to emphasize the regions within the disk photospheres.⁹ The figures also show that α has a strong correlation with the advective fraction, which increases as the surface density (and temperature) decrease on the upper branch. The reason for this trend of advective fraction is that much of the opacity on the upper branch is due to free–free absorption (i.e., Kramers opacity) and rises rapidly with falling temperature (∝T^−7/2), driving stronger convection.

Zoom In Zoom Out Reset image size

The reason why α is enhanced when thermal convection exists can be understood as follows. Convection drives vertical gas motions, which can create more vertical magnetic field than in the usual MRI turbulence. This assertion is supported by Figure 6(C), which shows that MRI turbulence is modified when convection exists so as to increase the time-averaged ratios of the vertical to the radial components for both velocity and magnetic field. Here the ratios are computed as

$$\begin{eqnarray} &f_{\rm mag} \equiv \frac{\int \left[\left<B_z^2\right>\right]\left[\left<p_{\rm thermal}\right>\right]dz}{\int \left[\left<B_x^2\right>\right]\left[\left<p_{\rm thermal}\right>\right]dz}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &f_{\rm vel} \equiv \frac{\int \left[\left<v_z^2\right>\right]\left[\left<p_{\rm thermal}\right>\right]dz}{\int \left[\left<v_x^2\right>\right]\left[\left<p_{\rm thermal}\right>\right]dz}. \end{eqnarray} \tag{ 20 }$$

The strengthened vertical magnetic field enhances the magnetic stresses since it is the seed for the most powerful, axisymmetric, modes of the MRI. Indeed, when convection acts, time-averaged and vertically integrated stress [W_xy] and the squared vertical magnetic field $[\int \langle B_z^2 \rangle dz]$ scale in proportion to one another, both deviating upwards from the trend ($\propto \bar{\Sigma }^{4/3}$) seen in the non-convective cases (Figure 6(D)). (The ratio of Reynolds stress to Maxwell stress hardly changes with $\bar{\Sigma }$ and is always ∼0.23.) In contrast, the integrated pressure [P_thermal] has a slightly decreasing trend when convection is strong. Apparently, the cooling due to convection lowers the time-averaged equilibrium pressure from the standard trend line despite the increased magnetic stress. Therefore, α, the ratio of stress to pressure, increases when convection is stronger.

In Figure 5, we also plot the profiles of the squared hydrodynamic Brunt–Väisälä frequency,

$$\begin{eqnarray} &\frac{N^2}{\Omega ^2} \equiv \frac{1}{\left[\left<\Gamma _1\right>\right]}\frac{d\ln \left[\left<p\right>\right]}{d\ln z} - \frac{d\ln \left[\left<\rho \right>\right]}{d\ln z}, \end{eqnarray} \tag{ 21 }$$

where Γ₁ ≡ (∂ln p/∂ln ρ)_s is the generalized adiabatic index and s is the specific entropy. Note that this hydrodynamic expression of buoyancy frequency is only valid where thermal pressure dominates magnetic pressure, i.e., the plasma β is larger than unity, which is interior to the two vertical dotted lines shown in the figure.

In solution (B), N² is negative in the midplane regions, precisely where the advective heat transport $\bar{F}^-_{\rm adv}(z)$ is substantial, indicating that the advective heat transport in the upper branch discussed in the previous subsection is really thermal convection associated with unstable hydrodynamical modes. We provide additional direct evidence that this is true thermal convection in the next subsection. Note that the large negative values of N² indicate that convection does not cause the time-averaged temperature gradient to be close to the adiabatic value, as supposedly happens in convection zones in stars. In fact, we have measured Mach numbers in the convective velocities as high as 0.1 in this simulation. Such superadiabatic gradients are also observed in MRI simulations with vertical convection by Bodo et al. (2013).

Apart from the case just mentioned, N² is generally positive, indicating convective stability, provided that β > 1. On the other hand, in solution (B), N² is negative in some regions where β < 1, which, however, does not mean that they are convectively unstable. Because magnetic pressure supports the plasma rather than thermal pressure alone, the hydrodynamic Brunt–Väisälä frequency is not the relevant quantity for buoyancy instabilities there. Rather, Parker instabilities can in general act in these regions (Blaes et al. 2007).

The two solutions near the right edge of the lower branch (ws0438F and ws0466F) also show a nonzero advective fraction (Figure 6(B)), but do not have an enhanced α (Figure 6(A)). We have confirmed that N² is negative where advective heat transport is observed in these simulations. However, the convection in these solutions is too weak to strengthen the vertical magnetic field, with convective Mach numbers that are two orders of magnitude smaller than on the upper branch.

So far we have been discussing trends in the time-averages of the simulations. However, in a given simulation, convection is often intermittent, and the system traverses limit-cycles, switching convection on and off episodically. These limit-cycles are traversed on roughly a thermal timescale.¹⁰ In contrast, the dwarf-nova limit cycles, because they require changes in the local surface density, are a separate phenomenon and take much longer, of order an inflow timescale.

An example of a convective/radiative limit cycle is seen in Figure 7, where time variations of the vertically integrated pressure $\tilde{P}_{\rm thermal}(t){ \equiv \lbrace P_{\rm thermal} \rbrace }$, stress $\tilde{W}_{xy}(t) {\equiv \lbrace W_{xy} \rbrace }$ and instantaneous advective fraction $\tilde{f}_{\rm adv}(t)$ defined as

$$\begin{eqnarray} &&\tilde{f}_{\rm adv}(t) \equiv \left\lbrace \frac{\int \left\lbrace \left<\left(e+E\right)v_z\right>\right\rbrace {\rm sgn}(z)\left<p_{\rm thermal}\right>dz}{\int \left\lbrace \left<\left(e+E\right)v_z+F_z\right>\right\rbrace {\rm sgn}(z)\left<p_{\rm thermal}\right>dz}\right\rbrace \end{eqnarray} \tag{ 22 }$$

are shown for Σ₀ = 140 g cm⁻² on the upper branch (ws0441F). Here the brackets { } denote box-car smoothing over a width of one orbit. Also shown are the instantaneous α and ratios of the vertical to the radial components for velocity and magnetic fields defined as

$$\begin{eqnarray} &{\tilde{\alpha }(t) \equiv \frac{\tilde{W}_{xy}(t)}{\tilde{P}_{\rm thermal}(t)}},\qquad\quad\ \ \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} &\tilde{f}_{\rm mag}(t) \equiv \left\lbrace \frac{\int \left<B_z^2\right>\left<p_{\rm thermal}\right>dz}{\int \left<B_x^2\right>\left<p_{\rm thermal}\right>dz}\right\rbrace, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &\tilde{f}_{\rm vel}(t) \equiv \left\lbrace \frac{\int \left<v_z^2\right>\left<p_{\rm thermal}\right>dz}{\int \left<v_x^2\right>\left<p_{\rm thermal}\right>dz}\right\rbrace. \end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

The curve of advective fraction demonstrates that convection occurs episodically, anti-correlated with the variation of pressure, indicating that convection is controlled by the temperature-sensitive opacity. The figure also shows that the ratios $\tilde{f}_{\rm mag}(t)$ and $\tilde{f}_{\rm vel}(t)$, whose time-averaged versions are enhanced when convection acts as discussed above, are actually enhanced at the beginning of each of the convective episodes. We interpret this as being due to the generation of vertical magnetic field by the onset of vertical convection, which seeds the axisymmetric MRI. The vertical to radial magnetic field ratio then falls back to the usual value as horizontal field is built by the MHD turbulence. The figure also shows that the stress begins to increase when the convection is fully developed and is followed by pressure with a finite time lag of several orbits. The stress parameter α, which is already higher than that of normal MRI turbulence, is further amplified when stress is high while pressure is low.

That thermal convection is actually operating during what we are calling the convective periods ($\tilde{f}_{\rm adv} \sim 1$) is clearly demonstrated by Figure 8, in which various quantities on a selected x–z plane (y = 0) at t = 103 orbits are shown. The specific entropy¹¹ is highest near the midplane, which drives low-density and high-temperature plumes that coherently transport heat upward (i.e., ev_z is mostly positive (negative) for z > 0 (z < 0), respectively). We therefore expect coherent vertical magnetic fields will be generated on the scale of the convective plumes, which is about half the pressure scale height.¹² Note that strong isolated magnetic fluxes are also associated with low density blobs, which suggests that the finite amplitude, slow mode buoyancy mechanism that we pointed out in Blaes et al. (2011) also contributes to vertical advection of heat. However, it is now completely dominated by the genuine thermal convection that fills much of the volume. These features in a convective period are contrasted with those in a radiative period at t = 90 orbits (Figure 9). The vertical gradient of specific entropy is now almost zero, indicating that the disk is convectively neutrally stable. This adiabatic gradient is caused by convection in the preceding convective period. Also, we see mostly random motions and only the slow mode mechanism is operating for the small net vertical advection of heat; however, the main heat transport mechanism here is, of course, radiative diffusion ($\tilde{f}_{\rm adv} \sim 0$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It might be instructive to visualize the convective/radiative limit cycles as trajectories in the pressure versus stress plane. Figure 10 shows such trajectories for Σ₀ = 140 (ws0441F), which we discussed above, as well as Σ₀ = 191 (ws0427F), 248 (ws0472F), 402 (ws0469F), and 1075 g cm⁻² (ws0429F) on the upper branch, in terms of the time variation of the vertically integrated pressure $\tilde{P}_{\rm thermal}(t)$ and stress $\tilde{W}_{xy}(t)$. The range of time is t₁ < t < t₂, where t₁ and t₂ are given in Table 1. The color intensity represents the instantaneous advective fraction $\tilde{f}_{\rm adv}(t)$. When the surface density is higher and convection is not present (Σ₀ = 402 and 1075 g cm⁻²), no limit cycle is seen, as pressure is almost unchanged in the face of stress fluctuations and thus the trajectories are almost vertical. As the surface density decreases (Σ₀ = 248, 191 and 140 g cm⁻²), the pressure fluctuation becomes larger and a limit cycle running clockwise in the plane is established: (1) magnetic stress is strengthened in fully developed convection, (2) stronger magnetic turbulence leads to greater dissipation, which increases the temperature and therefore the pressure, (3) higher temperature reduces the opacity, suppressing convection, and (4) without convection, magnetic fields weaken and the temperature declines, increasing the opacity, which eventually restores convection. On the left edge of the trajectories in the lower surface density (convective) cases, stress increases while pressure stays low; this phase lag further increases α as we discussed above.

Zoom In Zoom Out Reset image size

We always find runaway cooling (heating) of the disk beyond the left (right) edge of the upper (lower) branch, respectively. In Figure 11, we show time trajectories of such runs in the surface density versus effective temperature plane: ws0488R showing runaway cooling and ws0467R showing runaway heating. The initial development of MRI turbulence in both runs was similar to that in the other runs. However, both runs passed by the edge of the nearest stable branch and did not reach a steady state, which indicates that the two stable branches are actually truncated there. These facts indicate that any limit cycle in this plane runs anti-clockwise, which is consistent with the DIMs, where the cooling rate (heating rate) always exceeds the heating rate (cooling rate) beyond the left (right) edge of the upper (lower) branch, respectively. In fact, signs of the thermal runaway for the state transitions are seen near the edges of the branches. For example, the disk at the right edge of the lower branch (ws0466F) stayed in thermal equilibrium for ∼80 orbits and then began to flare up; on the other hand, the disk at the left edge of the upper branch (wt0442F) began to collapse after thermal equilibrium of ∼100 orbits. Similar behavior at the edges of the two stable branches was reported by Latter & Papaloizou (2012).

Zoom In Zoom Out Reset image size

Run ws0488R was stopped at t ∼ 50 orbits since the disk collapsed and resolving MRI near the disk midplane badly failed. Run ws0467R was also stopped at t  ∼  100 orbits when the mass loss from the simulation box became substantial (∼25%).¹³ The runs near the edges of the branches (ws0466F and wt0442F) were also stopped for the same reasons. To further simulate the thermal runaways or state transitions between the upper and lower branches, we would need to dynamically change the box size so that the MRI is always well resolved and the mass loss is kept small enough.

We could not find long-lived equilibria between the upper and lower branches, suggesting that the negative sloped portions of the alpha-based S-curves are unstable. In fact, in a few fiducial runs we found equilibria that lasted more than several tens of orbits, but these were not fully reproduced when the box size or the resolution was changed.

To check the robustness of our results, we performed two kinds of tests; one to check dependence on the initial conditions (Figure 12) and the other to check numerical convergence (Figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For four arbitrarily selected fiducial runs, we ran a supplementary simulation whose surface density is almost the same, but whose initial effective temperature is different from the fiducial run. In Figure 12, the initial and final states of the four pairs of fiducial and supplementary runs (ws0441F and ws0494C, ws0464F and ws0465C, ws0471F and ws0472C, and ws0468F and ws0469C) are shown. The final states of paired runs are roughly the same both in effective temperature and α, which confirms that the two branches are actually unique attractors.

We also arbitrarily selected five fiducial runs on the upper branch, Σ₀ = 1075 (ws0429F), 402 (ws0469F), 247 (ws0471F), 179 (ws0437F), and 132 g cm⁻² (ws0446F), and for each we ran two supplementary simulations: doubling the horizontal box size in one, and increasing the resolution by 1.5 times and the box size by 1.2 times in the other one. As shown in Figure 13, the final states of the fiducial run and the corresponding two supplementary runs are roughly the same, both in effective temperature and α for all five cases. We may therefore conclude that our two main results, the thermal equilibrium curve and the high α near the edge of the upper branch, are not sensitive to the box size or the resolution.

As shown in Figure 3, those simulations that do not exhibit convection lie on our DIM-model curve with fixed α = 0.03 on the upper branch, while they are slightly below the curve on the lower branch (cf. Latter & Papaloizou 2012). We suspect that the discrepancy may come from our DIM's assumption that the disk is very optically thick even on the lower branch, where the measured optical depth in the simulations can be as low as 2.3 (Table 1).

On the other hand, our results near the left edge of the upper branch deviate upward from the fixed α curve as α increases due to convection. In some cases, however, our results with larger values of α lie below a DIM-model curve computed assuming a smaller α. Also, the minimum (maximum) surface density of our upper (lower) branch, respectively, are larger than those of the relevant DIM-model curve. These discrepancies are presumably due to our neglect of convection in computing the DIM curves because convection tends to increase the critical surface densities (see, for example, Pojmański 1986).

It is widely known that the saturation of MRI turbulence in the local shearing box depends on the net vertical flux or net toroidal flux threading it (see, for example, Hawley et al. 1995; Latter & Papaloizou 2012). Although we do not impose a net vertical flux in the box, we do impose a net toroidal flux initially in the box (see Section 2.4). Therefore one might argue that the saturation of MRI turbulence in our simulations could be affected by the initial net toroidal flux.

Stratified shearing box simulations, however, generally do not retain a net toroidal flux because of buoyant escape through the vertical boundaries. Figure 14 shows the time variation of the net toroidal flux $\tilde{\Phi }_y(t)\equiv \lbrace \int \langle B_y \rangle dz \rbrace$ as well as the net radial flux $\tilde{\Phi }_x(t)\equiv \lbrace \int \langle B_x\rangle dz \rbrace$ for the initial hundred orbits for three selected runs: a convective solution on the upper branch (ws0441: α = 0.0927), a radiative solution on the upper branch (ws0429: α = 0.0332), and a radiative solution on the lower branch (ws0435: α = 0.0312). The net toroidal flux in every case fluctuates significantly in time, flipping in sign, and there is no indication of any memory of the initial net toroidal flux. (The reason for the sign flips is presumably that azimuthal flux arises from shear of radial flux, which also flips sign over time due to the still poorly understood dynamo of stratified MRI turbulence.) Therefore, we conclude that the saturation of MRI turbulence in our simulations is independent of the initially imposed toroidal flux.

Zoom In Zoom Out Reset image size

We have shown that convection enhances MRI turbulent stress, which can increase α above 0.1. Since convection necessarily occurs in the outburst phase due to the strong temperature dependence of opacity, we have found an α-enhancement mechanism that is due to the internal physics within the disk in this regime.

There are other external mechanisms, however, that might be considered candidates for producing enhanced α. For example, it is well known that imposition of net vertical magnetic flux raises the saturation level of the turbulence. However, it is also known that the dependency of stress on net vertical field B_znet is fairly strong, ∝B_znet² (Suzuki et al. 2010; Okuzumi & Hirose 2011). Therefore, the fact that the observed α in the outburst phase is always of order 0.1 is a strong constraint on the existence of a global net vertical field in the disk in the outburst phase (King et al. 2007). On the other hand, global simulations of MRI turbulent disks have produced local net vertical fluxes through magnetic linkages in the disk corona (Sorathia et al. 2010). It could be that this mechanism might play a role in producing large α's in the outburst phase, but why this would only occur in the outburst phase is unclear. Sorathia et al. (2012) have also suggested that the large α values inferred in the outburst phase may be due to transient periods of magnetic field growth in the jump to outburst, together with gradients in the global disk. For these and many other reasons, thermodynamically consistent, global MRI simulations of disks in the hydrogen ionizing regime will be of great interest.

Another point worth mentioning is that an enhanced stress does not necessarily lead to an enhanced α. If, for example, pressure also rises in proportion to the enhanced stress (via enhanced dissipation), α, the ratio of stress to pressure, would not be increased. In our simulations, the convective cooling controls the pressure well enough to lead to an increase in α. Net vertical flux, which increases stress, could in principle also explain the high α in the high state, but whether or not it does will depend on the scaling of cooling rate with pressure in the presence of vertical flux. New simulations that carefully account for thermodynamics will be necessary to determine this scaling.

We have remarked on the possible thermal instability of the S-curve branch that should link the low and high states. Any such instability would be qualitatively different from any thermal instability that affects a radiation pressure dominated regime (Shakura & Sunyaev 1976; Turner 2004; Hirose et al. 2009; Jiang et al. 2013). In the temperature range between the low and high states relevant to dwarf novae, the opacity increases rapidly with increasing temperature because this is the regime in which H ionizes, so that the dominant opacity source is free–free or bound–free absorption of the negative hydrogen ion H⁻. Thus, an upward fluctuation in heating receives positive feedback because it is accompanied by weaker cooling. By contrast, thermal fluctuations in the radiation-dominated regime are aided both by the sensitivity of radiation pressure to temperature (∝T⁴) and possible dynamical coupling between total thermal pressure (gas plus radiation) and heating associated with MHD turbulence (Jiang et al. 2013).