{'fraction_non_alphanumeric': 0.02971809769743921, 'fraction_numerical': 0.01719388853023456, 'mean_word_length': 4.50539035659282, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Thermal Stability of Radiation Dominated Accretion Disks

We study the long-term thermal stability of radiation dominated disks in which the vertical structure is determined self-consistently by the balance of heating due to dissipation of MHD turbulence driven by the magneto-rotational instability (MRI), and cooling due to radiation emitted at the photosphere. The calculations adopt the local shearing box approximation, and utilize the recently developed radiation transfer module in the Athena MHD code based on a variable Eddington tensor rather than an assumed local closure. After saturation of the MRI, in many cases the disk maintains a steady vertical structure for many thermal times. However, in every case in which the box size in the horizontal directions is at least one pressure scale height, fluctuations associated with MRI turbulence and dynamo action in the disk eventually trigger a thermal runaway which causes the disk to either expand or contract until the calculation must be terminated. During runaway, the dependence of the heating and cooling rates on total pressure satisfy the simplest criterion for classical thermal instability. We identify several physical reasons why the thermal runaway observed in our simulations differ from the standard alpha disk model, for example the advection of radiation contributes a non-negligible fraction to the vertical energy flux at the largest radiation pressure, most of the dissipation does not happen in the disk mid-plane, and the change of dissipation scale height with mid-plane pressure is slower than the change of density scale height. We discuss how and why our results differ from those published previously. Such thermal runaway behavior might have important implications for interpreting temporal variability in observed systems, but fully global simulations are required to study the saturated state before detailed predictions can be made.


INTRODUCTION
The inner regions of black hole (BH) accretion disks are thought to be radiation pressure dominated whenever the accretion rate is larger than a few percent of the Eddington limit. Thus, understanding the properties of radiation dominated accretion flows is essential in order to be able to interpret and predict the spectrum from luminous systems such as quasars. Shortly after the standard α disk model was proposed (Shakura and Sunyaev 1973), it was realized that radiation dominated accretion disks are thermally unstable if the stress is proportional to the total pressure P t (Shakura and Sunyaev 1976). Piran (1978) proposed a general criterion to determine thermal stability based on the dependence of the total cooling (Q − ) and heating (Q + ) rate on P t : the disk is thermally unstable if: where Σ is the surface density. In a radiation dominated α disk model with electron scattering as the dominant source of opacity and stress τ rφ = αP t , then Q + ∝ P 2 t while Q − ∝ P t at fixed surface density (Piran 1978), thus the inequality (1) is satisfied and the disk is thermally unstable. The stability of radiation-pressure dominated disks is quite different from thermal instability in the much cooler disks associated with dwarf novae and low-mass X-ray binaries (Lasota 2001). The former instability is driven by strong variations in heating, while the latter by strong variations in the radiative cooling due to the very sharp opacity shifts near 5000K. There are, however, several uncertainties in the α disk model which cast into doubt whether the thermal instability actually exists in real disks. First and foremost is the assumption that the stress is proportional to total pressure. For example, if instead the stress is assumed to be proportional to gas pressure alone (the so-called β model), the disk is thermally stable (Sakimoto and Coroniti 1981;Stella and Rosner 1984;Merloni 2003). More importantly, it is now realized that the magneto-rotational instability (MRI, Balbus and Hawley 1991) is the physical mechanism which generates the stress in the inner (fully ionized) regions of black hole accretion disks (Balbus and Hawley 1991). Dissipation of the turbulence produced by the MRI results in heating. However, there is no reason to expect the heating rate produced by the MRI should be proportional to the total pressure, nor that the vertical profiles of dissipation and heating should be the same as in an α disk model (e.g., Zhu and Narayan 2013). Thus, the thermal stability of MRI unstable disks is uncertain, and can only be investigated with nonlinear radiation MHD simulations which capture the MHD turbulence, heating and cooling of the disk self-consistently. Evidence of dwarf novae thermal instability and associated limit cycle behavior have been found from MRI simulations (Latter and Papaloizou 2012).
Recently, the thermal properties of radiation dominated MRI unstable disks have been investigated in a series of papers (e.g., Hirose et al. 2009), using a module for radiation transport based on flux-limited diffusion (FLD) developed for the ZEUS MHD code (Stone and Norman 1992) by Turner and Stone (2001). Using 3D local shearing-box simulations of stratified MRI unstable disks that extend for as long as ∼ 60 thermal times, they report the important result that radiation pressure dominated disks are thermally stable (Hirose et al. 2009). Moreover, they develop a simple model to interpret this result, based on the observed time lag between stress and pressure fluctuations which indicates that stress is not determined by total pressure.
In this paper, we report results from a new investigation of the MHD of radiation dominated disks using more advanced numerical methods ). While we have been able to reproduce many of the pioneering results reported by Hirose et al. (2009), in one important respect our results differ: in our calculations the disk always eventually undergoes thermal runaway (either expansion or collapse). Working directly with Hirose et al., we are only able to reproduce their results if we also adopt FLD and use the same horizontal box size as well as the initial condition as used in Hirose et al. (2009). Whenever the horizontal box size is at least one scale height wide, or the more accurate radiation transfer algorithm is used, the thermal runaway shows up. The similarities and differences between the results from the two codes will be reported in a future joint publication (Jiang et al., in preparation). In this paper, we simply report the thermal evolution of radiation pressure dominated disks observed in our new simulations. This paper is organized as follows. The equations we solve are given in § 2, and the initial and boundary conditions we use in § 3. Our primary results are described in § 4, while § 5 presents a discussion and conclusion.

EQUATIONS
We solve the equations of radiation MHD in a frame rotating with orbital frequency Ω at a fiducial radius r 0 from the central BH, using the local shearing box approximation (Hawley et al. 1995;Hirose et al. 2009). Curvature of the orbit is neglected so the equations are solved in the local Cartesian coordinate (x, y, z) with unit vectors (î,ĵ,k), which represent the radial, azimuthal and vertical directions respectively. With the vertical component of the gravitational force of the central black hole and the energy exchange due to Compton scattering (Hirose et al. 2009) added, the equations are (e.g., ) where the radiation source terms are, In the above equations, the shear parameter q ≡ −d ln Ω/d ln r is 3/2 for a Keplerian disk. The pressure P * ≡ (P + B 2 /2)I (with I the unit tensor), and the magnetic permeability µ = 1. The total gas energy density is E = E g + ρv 2 /2 + B 2 /2, where E g = P/(γ − 1) is the internal gas energy density with γ = 5/3. The radiation constant a r = 7.57 × 10 15 erg cm −3 K −4 . The radiation temperature T r is defined as T r ≡ (E r /a r ) 1/4 , while T e is the temperature equivalent to the electron rest mass, T e = 5.94 × 10 9 K. The frequency mean absorption and scattering opacities (attenuation coefficients with unit cm −1 ) are denoted by σ aF and σ sF , while σ aP and σ aE are the Planck and energy mean absorption opacities.
Equations 2 and 3 are solved in the mixed frame, that means the radiation flux F r and energy density E r are Eulerian variables, while the material-radiation interaction terms in equations 3 are written in the comoving frame (e.g., Lowrie et al. 1999). The Eulerian and co-moving flux F r,0 are related through F r,0 = F r − (vE r + v · P r ). The radiation pressure P r and energy density are related through a variable Eddington tensor (VET), P r = fE r . We calculate f from a formal solution of the transfer equation using the method of short characteristics . For the calculations presented in this paper, we use 10 angles per octant to compute the VET.
The equations are solved with the new Godunov radiation MHD code as described and tested in  and , with the improvements given by Jiang et al. (2013). An orbital advection scheme (Stone and Gardiner 2010) is used to speed up the calculation. The Compton heating term is separated from other terms and added to the energy equation in the same way as done by Hirose et al. (2009). These ideal MHD equations do not include explicit viscosity or resistivity. Instead, dissipation occurs at the grid scale through numerical diffusion. However, because we solve the total energy equation with a conservative method, the kinetic and magnetic energy dissipated at the grid scale is transformed into thermal energy, and is not lost.

INITIAL AND BOUNDARY CONDITIONS
In order to facilitate comparisons between our results, we adopt the same values for the physical parameters as Hirose et al. (2009), as given in Table 1 of their paper. 1 The local frame of reference represented by the shearing box is centered at r 0 = 30 GM/c 2 , with the mass of the central black hole M = 6.62M ⊙ , and corresponding orbital frequency Ω = 190.1 s −1 . The mean molecular weight is assumed to be 0.61. Following Hirose et al. (2009), we include electron scattering opacity 0.33 cm 2 g −1 , Plank-mean free-free absorption opacity 3.7 × 10 53 ρ 9 /E 7 g 1/2 cm 2 g −1 and Rosselandmean free-free absorption opacity 1.0 × 10 52 ρ 9 /E 7 g 1/2 cm 2 g −1 .
The initial vertical profile of the disk is calculated in the same way as described in Section 2.4 of Hirose et al. (2009). For the magnetic field, we initialize two oppositely twisted flux tubes with the same net azimuthal flux within −0.25 < z < 0.25. The B x and B z components of the field are generated by the vector po- z 1/2 for |z| < 0.8. This field configuration is slightly different from that used by Hirose et al. (2009) (it is symmetric with respect to the disk midplane), however we have confirmed that different choices for the initial magnetic field do not change our results on the thermal stability of the disk. But different initial vertical profiles can affect how fast the thermal runaway happens. We apply a density floor 5 × 10 −6 of the initial mid-plane density throughout the numerical integration to avoid small time steps when the density in some cells become too low.
Shearing periodic boundary conditions (Jiang et al. 2013) are used for the x-direction, and simple periodic for the y-direction. For the vertical direction, the gas pressure and density are copied to the ghost zones from the last active cells. The three components of the velocity are also copied when the vertical component is outward. To prevent inflow, they are set to zero when the vertical component is inward. For the magnetic field, we copy all three components to the ghost cells when the velocity is outward, and set the horizontal components to zero (copying only the vertical component) when the velocity is inward. We have also tried to always copy the three components of magnetic field to the ghost cells even velocity is inward, with no effect on the results presented here. Anomalous resistivity (e.g., Sano et al. 2004) is used within one scale height from the vertical boundary to overcome the numerical difficulties in the strongly magnetized photosphere as described in Hirose et al. (2009). For the radiation field, we assume zero incoming specific intensity in the short characteristics solver used to compute the VET. To be consistent, we set |F r,z |/ (cE r ) = √ f zz in the ghost zones, where f zz is the zz component of the Eddington tensor in the ghost zones, that is we assumes strictly outgoing radiation flux. The horizontal components of the radiation flux are copied to the ghost zones from the last active zones. The radiation energy density E r in the ghost zones is then calculated from ∂(f zz E r )/∂z = −(σ sF + σ aF )F r,z .

RESULTS
From the laminar initial conditions, vigorous MHD turbulence is produced by the MRI within a few orbits. Thus, the central density and temperature, as well as the vertical profiles of these quantities, quickly evolve away from their initial values. After saturation of the MRI, the disk is heated by the turbulence and cooled by radiation from the photosphere. As reported by Hirose et al. (2009), we find that after saturation of the MRI the disk can maintain an approximate equilibrium for many thermal times. However, we also find that for all the numerically converged simulations, inevitably the disk undergoes a thermal runaway, either heating up and expanding until the photosphere hits the boundary of our domain, or cooling down and collapsing until it can no longer be resolved by our numerical mesh. We investigate this behavior with a series of calculations using different surface densities described below. Parameters and outcome of the simulations are summarized in Table 1.

Case A: Large Surface Density
We first study a case with large surface density (twice the value used by Hirose et al. 2009), giving an electron scattering optical depth across the whole disk τ es = 7.11 × 10 4 . The parameters for this run, hereafter referred to as RHVET, are listed in Table 1. The fiducial units used in this calculation are ρ 0 = 0.0566 g cm −3 , T 0 = 2.45 × 10 7 K, and P 0 = 1.89 × 10 14 dyn cm −2 . The ratio of radiation-to-gas pressure at the mid-plane initially is 4.13.
We first ran the calculation adopting the Eddington approximation, f = 1/3I, to save computer time. The MRI saturates within the first 10 orbits, and the disk maintained a steady structure until roughly 50 orbits. The total (Maxwell plus Reynolds) stress in the saturated state normalized by the total pressure is 0.023. Thereafter, we observed the disk to undergo a runaway, heating up and expanding until the total pressure increased by a factor of ∼ 6. Within the first 150 orbits of evolution, only 1.5% of the total mass is lost through the open vertical boundaries.
Concerned that this behavior is due to the approximate treatment of radiative transfer, we restart the simulation at 50 orbits using a VET computed with short characteristics. Figure 1 shows space-time plots of the vertical profiles of the horizontally averaged density and azimuthal magnetic field from this calculation. The position of the photosphere for electron scattering opacity is shown as the white line in the top panel. It corresponds roughly to the point at which the density drops to 10 −4 ρ 0 . Note the photosphere is well inside the simulation domain for t < 100 orbits, but has reached the top and bottom of the box when we terminate the simulation. The thermal time, defined as the ratio between the total energy within the simulation domain and the total cooling rate, gradually increases from 10 to 20 orbits. Figure 1 shows that after an initial period of 50 orbits, once again the disk undergoes a runaway over the subsequent ∼ 10 thermal times.
The space-time plot of the toroidal magnetic field B y shows reversals roughly every 10 orbits, reproducing the well-known butterfly diagram for the MRI in stratified disks, which is observed both with Hirose et al. 2009;Blaes et al. 2011) and without (Stone et al. 1996;Davis et al. 2010) radiation. This pattern is driven by a dynamo process in the disk (Brandenburg et al. 1995;Gressel 2010;Blackman 2012). The photosphere of the disk also moves up and down in concert with this pattern. Note the toroidal magnetic field fluctuations become noticeably stronger as the disk gets thicker towards the end of the simulation. The volume averaged Maxwell and Reynolds stress also increase secularly as the disk is heated up, increasing by a factor of ∼ 3 by the end of the simulation. This increase is consistent with the scaling of the stress due to the MRI with domain size first noticed by Hawley et al. (1995); the larger the disk thickness, the larger the size of turbulent eddies driven by the MRI, and the larger the stress.
To test whether the behavior we observe is consistent with the basic tenets of thermal instability, we calculate the heating and cooling rates as follows. The heating rate per unit area Q + in a shearing box simulation of the MRI is set by the work done on the fluid by the radial boundaries (Hawley et al. 1995;Gardiner and Stone 2005), therefore where the integration is one radial side of the domain, and dy and dz are the cell size along the y and z directions respectively. The cooling rate per unit area Q − is dominated by the radiation flux leaving the top and bottom of the simulation domain, therefore We plot Q + and Q − versus the mid-plane pressure P z0 , as well as the change of P z0 with time, in Figure 2. Consistent with the fact that the mid-plane pressure always increases with time, Q + is always larger than Q − . De-spite large amplitude fluctuations in Q + and Q − due to the chaotic nature of MRI turbulence, the mean values show a clear trend, increasing with P z0 . We fit a powerlaw to both over the period 80 − 130 orbits. During this time, Q + and Q − start from almost the same value, the photosphere is still well inside the simulation domain, and the space time plots clearly show the thermal runaway occurring. The best fit gives Q + ∝ P 1.60 z0 while Q − ∝ P 0.98 z0 during this time. Thus, our result is consistent with the most basic premise of thermal instability (equation 1), namely the change of the heating rate with pressure is faster than the change of cooling rate. However, the difference between the power law indices of Q + and Q − is smaller than the predicted value of traditional linear instability as discussed in the introduction.

Case B: Moderate Surface Density
Next, we describe the results of another simulation which uses a surface density identical to that used by Hirose et al. (2009), that is one half of the value used in RHVET described above, giving an electron scattering optical depth across the whole disk τ es = 3.56 × 10 4 . Hereafter we refer to this simulation as RMLVET. Parameters for this calculation are given in Table 1. The fiducial units used in the calculation are ρ 0 = 0.0566 g cm −3 , T 0 = 2.45 × 10 7 K, and P 0 = 1.89 × 10 14 dyn cm −2 . The ratio of radiation-to-gas pressure at the midplane initially is 4.80. Note that while the parameters and initial conditions for this calculation are identical to Hirose et al. (2009), the computation is performed in a domain which is more than two times larger in each dimension.
Once again, we first evolve the disk adopting the Eddington approximation f = 1/3I for the first 20 orbits and restart the simulation with VET. After saturation of the MRI, and for the first 60 orbits of this simulation, the total energy density and Maxwell stress fluctuate around mean values with no systematic trend, similar to the results reported by Hirose et al. (2009). The ratio between the total stress and pressure during this time is 0.034. However, once again after 60 orbits the energy density and stress being to increase systematically, and the disk expands, albeit with large fluctuations in both during this period. During the expansion phase, our best fit power laws to the heating and cooling rates give Q + ∝ P 2.42 z0 while Q − ∝ P 1.44 z0 , which is again consistent with the thermal instability criterion (equation 1).
Interestingly, if we continue the simulation with the Eddington approximation, we find the disk collapses after 80 orbits instead of expanding.

Case C: Low Surface Density
Finally, we report the evolution of a disk with low surface density, one quarter the value used by Hirose et al. (2009) giving a total optical depth τ es = 8.36 × 10 3 . Parameters for this calculation, hereafter referred to as RSVET, are given in Table 1. Initial parameters of this simulation are chosen according to the radiation pressure dominated solution of the α disk model with assumed α = 0.0125 and 42.5% of the Eddington accretion rate. using fiducial units ρ 0 = 0.00212 g cm −3 , T 0 = 2.87 × 10 7 K, and P 0 = 8.29 × 10 12 dyn cm −2 . The ratio of radiation-to-gas pressure at the mid-plane initially is 206. We use our VET module to compute the entire evolution. The space-time diagram for this run after the initial 10 orbits is shown in Figure 3. For this simulation, the stress increases for the first 20 orbits, during the saturation phase of the MRI, reaching a value of 0.017 of the total pressure. Thereafter, the stress, radiation, and gas internal energy densities all start to decrease. Unlike the previous cases, this strongly radiation pressure dominated disk does not show any period of thermally stable structure. After another 20 orbits of evolution, the disk has collapsed to a degree that we can no longer resolve the MRI. During the collapse, the best fit power laws to the heating and cooling rates give Q + ∝ P 1.90 z0 while Q − ∝ P 0.90 z0 . Therefore this solution is also consistent with the thermal instability criterion in that when the pressure drops, the heating rate decreases faster than the cooling rate.

The Effects of the Radiation Transport Algorithm
and Domain Size In order to assess the impact of the radiation transfer algorithm on our results, we have implemented a FLD module in Athena. 2 We first consider the case of low surface density and high radiation to gas pressure ratio. We restarted from a ZEUS simulation with the same parameters as the RSVET simulation reported in Table 1 (kindly provided by S. Hirose), which is labeled RSFLD. The initial cooling rate with FLD in Athena is the same as in the ZEUS solution. The disk maintains a steady structure for ∼ 20 orbits. However, after 20 orbits the solution still collapses in a similar manner to original RSVET simulation reported above. If we restart the same ZEUS solution with our VET module or Eddington approximation, the initial cooling rate is increased by a factor of ∼ 2, indicating that FLD underestimates the cooling rate relative to VET. As a result, the disk cannot maintain a steady structure and collapse proceeds immediately. These results suggest that the Athena simulations always undergo a thermal runaway, regardless of the radiation transfer algorithm, for these low surface density runs where the ratio of radiation-to-gas pressure is larger than ∼ 200.
In domains with larger surface density, the radiationto-gas pressure ratio is smaller (∼ 10 for RMLVET). In 2 We note that the Athena and ZEUS implementations of FLD are not identical, primarily owing to the differences in the underlying methods (finite volume versus operator splitting) and the choice of frame used to specify the radiation transfer equations and variables (co-moving versus mixed-frame). Heating (red line) and cooling (black line) rates versus mi-dplane pressure for RHVET. Also plotted is the evolution of the mid-plane pressure with time (blue line). The red and black dashed lines are the best fit power-law between 80 − 130 orbits. The unit for Pt is P 0 , while the units for Q + and Q − are 0.21H 3 P 0 Ω. order to assess the impact of the radiation transfer algorithm and the impact of domain size in this regime, we run two simulations labeled RMEdd and RMFLD, utilizing the Eddington approximation and the FLD algorithm, respectively. These use the same simulation parameters as RMLVET except that both L x and L y are smaller by a factor of 2, to match the 1112a simulation of Hirose et al. (2009). We find that the disk collapses after ∼ 150 orbits in RMEdd. In contrast, RMFLD shows no clear thermal runaway over 350 orbits, consistent with the behavior described in Hirose et al. (2009). While the disk shows periods of expansion, it eventually falls back to its original structure. The heating and cooling rate track each other very well during the simulation. Therefore, Athena and ZEUS agree reasonably well for these higher surface density runs, when they utilize the same horizontal domain size and radiation transfer algorithm (FLD). However, when the more accurate radiation transfer algorithm is used, the disk is still undergoes a thermal runaway even with a smaller horizontal box size.
To further examine the effect of the horizontal box size for the case when radiation pressure is only 10 times the gas pressure, we consider another run labeled RMFLDL, where we double the horizontal box size of RMFLD run so that L x = H and L y = 4H while keeping all other parameters to be the same. This run has the same simulation parameters as RMLEdd and RMLVET but adopts FLD. The disk collapses within 50 orbits in this case, which is consistent with the thermal runaways observed in RMLEdd and RMLVET, but different from RMFLD. This experiment shows that by just increasing the hori- zontal box size from ∼ 0.5H × 2H to H × 4H, the disk can change from the stable behavior to be unstable, even when FLD is used.
The effects of the horizontal box size on the saturation state of MRI have been studied in previous isothermal MRI simulations. Simon et al. (2012) studied the saturation of MRI through a series of simulations with different horizontal box sizes. They found that the fluctuation of MRI turbulence decreased with increasing horizontal box size. Furthermore, properties of the MRI turbulence, such as vertical dissipation profile and correlation length, are only converged when the horizontal box size is larger than H × 4H. It has also been observed that small box MRI simulations shows large spikes in stress, and thus heating rate, due to the recurrent channel solutions, which are significantly reduced when the horizontal box size is increased (Bodo et al. 2008;Jiang et al. 2013). Larger box size also increases the number of modes in the system and phase trajectories of the system are less constrained. These results suggest that the simulations with larger domains are more reliable to determine whether the radiation dominated disks are thermally stable or not.
These results indicate the outcome of the simulations depends on the radiation-to-gas pressure ratio, the box size, and the radiation transport algorithm used. For modest radiation-to-gas pressure ratios (∼ 10), the box size and the algorithm used both have an impact, although our preferred setup (larger domains with the VET algorithm) produces thermal runaway. When the radiation pressure is ∼ 200 times the gas pressure and the surface density is small, the differences caused by the box size and radiation transfer algorithm are less important and we observe a thermal runaway for all the simulations.

DISCUSSION
The final outcome of all the simulations reported above is summarized in Table 1. We will report a detailed analysis of these simulations, as well as the results from many more that explore a much wider range of radiation-togas pressures, in a future paper (Jiang et al., in prep.). Here, we first show that although we find thermal runaway as classical thermal instability predicted, many assumptions in the α disk model are not satisfied in our simulations and the details of the thermal runaway differ from the predictions.Then we discuss how our radiation dominated results compare to simulations run with the ZEUS code.

Comparison with the α-disk Model
Before we begin an in-depth comparison, it is useful to state precisely what is meant by the α-disk model in the current context. In terms of vertical structure, the αdisk is effectively a "one-zone" model as all variables are characterized by a single value at a given radius. Since we are comparing with local simulations with fixed mass, it is useful to consider the properties of the disk at a single radius, with constant values of surface density Σ, Ω, and surface radiation flux F rs . One to one correspondence is assumed between stress and the total mid-plane pressure τ rφ = αP z0 , with α constant. We focus here on the radiation dominated limit where P r (0) ≡ P z0 .
The disk is assumed to be in hydrostatic and thermal equilibrium, which given Σ, Ω, and F rs determines the mid-plane radiation pressure via Here κ s is the electron scattering opacity and ξ is a parameter that depends on the vertical distribution of dissipation, which is assumed to be constant in the α-disk model. If one follows (Shakura and Sunyaev 1973) and assumes dF r /dz ∝ ρ, ξ = 1/2. The cooling rate per unit area in the disk is assumed to be determined solely by radiative cooling so that Q − = F rs . The heating rate per unit area is proportional to an integral over z of the stress with H s defining a characteristic scale height for the stress. Finally, hydrostatic equilibrium determines a characteristic flux scale height Thermal instability follows from equation 1 by assuming ξ and α are constants and that H s = H F . We then find that Q + ∝ P 2 z0 and Q − ∝ P z0 . This simple model then predicts a linear instability that grows on the order of the thermal time ≃ 1/(Ωα).
Strictly speaking, none of these assumptions are obeyed by the shearing box simulations, therefore it is not surprising that we do not find an exponential runaway on the thermal timescale. For the moderately radiation pressure dominated case, the disk does in fact stay in a roughly equilibrium state with large amplitude fluctuations for several thermal times. Both RHVET and RMLVET survive for ∼ 50 − 60 orbits before thermal runaway takes over. In fact, it is not even clear how a classical linear analysis can be applied to an MRI turbulent disks, because the disk already contains nonlinear amplitude fluctuations in the stress, heating and cooling rates driven by the turbulence. Even the assumption of a single equilibrium state is questionable because variations in the vertical distributions of dissipation and stresses can, in principle, offer a range of equilibrium (or near-equilibrium) configurations to the turbulent flow.
Instead, the behavior we observe is perhaps more akin to a nonlinear instability, in which a finite amplitude fluctuation is required to drive the system away from equilibrium. In this interpretation, the equilibrium state in the radiation pressure dominated case may be better described as a basin of attraction which temporarily traps solutions. Even when there is no runaway, turbulence provides large amplitude fluctuations in the disk structure which causes the system to wander within the basin (Covas et al. 1997;Ashwin and Rucklidge 1998;Janiuk and Misra 2012). Exploring the topology of solutions in phase space for this dynamical system to discover what conditions fluctuations must satisfy to produce runaway is well beyond the goals of this paper. Nonetheless, such an interpretation may provide an explanation for the effects of horizontal box size on the thermal instability we find.
We can quantify the difference with the α-disk by computing some appropriately averaged quantities. We start by computing the scaleheight H s where all integrals are performed over the full extent of the domain. We plot H s and H F in Figure 4 for the simulations RHVET and RMEdd versus the total mid-plane pressure. We see that H F is not directly proportional P z0 due to fluctuations. It is also generally larger than the α model prediction (especially in simulation RMEdd) because the dissipation rate per unit mass from MRI turbulence actually increases with decreasing column density. This leads to more dissipation near the surface and a lower ratio of surface flux to mid-plane pressure. Especially when the density profile changes significantly with mid-plane pressure, H F deviates from linear dependence on P z0 significantly. The red lines shows best-fit power-law relation between the scale heights and pressure, fit from 80-130 orbits for RHVET and 20-150 orbits for RMEdd. In the RHVET run, this is close to linear for H F , suggesting that the profile of the dissipation per unit mass is relatively constant (at least when averaged over fluctuations), in agreement with the α-disk assumption. However, simulation RMEdd shows a tendency for a steeper than linear increase in H F with P z0 , indicating a a larger share of the dissipation occurs at low densities as the disk expands. We also see that H s differs significantly from H F , particularly at large pressure, due to a much weaker than linear increase with pressure (H s ∝ P 0.36 z0 for RHVET and H s ∝ P 0.66 z0 for RMEdd). This is a fairly general result that holds in almost all of our simulations.
The different dependence of H F and H s on pressure is closely related to the presence of vertical fluxes of energy other than radiative diffusion. During thermal runaway, the advective flux of radiation energy near the mid-plane due to turbulence and magnetic buoyancy (Blaes et al. 2011) is not negligible, as is assumed in the α disk model. Figure 5 shows vertical profiles of the temporally and horizontally averaged diffusive, advective, and Poynting flux for RHVET at two different times. Both the advective and Poynting flux peak at ∼ 2H. Thereafter both drop quickly to zero approaching the photosphere. The ratio between the advective and diffusion flux near the mid-plane becomes larger when the disk becomes hotter. At 90 orbits, the advective flux is larger than the diffusive flux within ∼ 2H. The diffusive flux is constrained by the vertical component of gravitational acceleration, which does not change, whereas the advective flux depends on the amplitude of the turbulence and magnetic energy (which determines buoyancy), both of which are increasing as the disk gets thicker. Note that at the photosphere, radiation is indeed the dominant cooling mechanism of the disk, as assumed in the α model: at the top and bottom of the simulation domain, the mechanical and Poynting fluxes are only ∼ 1% of the radiation flux in RHVET. Because of the extra cooling provided by the advective flux at the mid-plane, the radiative flux, and therefore, H F can increase without a proportional rise in H s .
The estimation of α is somewhat more arbitrary. We could, in principle, compute α as the ratio of the volume averaged stress to the volume averaged total pressure, but this would be contrary to the α model assumption, which implicitly scales α with the mid-plane pressure. The distinction is relevant because the vertically averaged pressure generally scales as the square of the midplane pressure. Simply using the ratio of the mid-plane stress and pressure is also problematic because the stress and pressure show markedly different profiles, even when time averaged. While the pressure peaks at or near the mid-plane, the stress peaks somewhat off the mid-plane with a local minimum at the mid-plane. Hence we prefer to define α = τ rφ /P z0 , where τ rφ is the volume average of the stress for |z| < H s .
For this definition, we find that α is generally a weakly increasing function of mid-plane pressure, but one with large fluctuations about the best-fit mean relations. For example, the RHVET simulation provides a best-fit α ∝ P 0.19 t0 . The this trend of weakly increasing α with pressure is a fairly general result in our simulations, and is in approximate agreement with the α model assumption. Along with the weak increase of H s with pressure, this implies a value of ∂ ln Q + /∂ ln P z0 < 2 for most simulations. In the specific case of RMFLD simulation, which shows no evidence for runaway, the dependence is sufficiently weak that ∂ ln(Q + /Q − )/∂ ln P z0 is only about 0.1.
Despite the differences discussed above, we have shown that during runaway our solutions agree with the most basic criterion for thermal instability: the heating rate Q + has a steeper dependence on mid-plane pressure P z0 than does the cooling rate Q − , albeit with significant fluctuations on shorter timescales. However, the actual scalings we measure can be substantially different from what the α model predicts. In particular, the differ-ence in the logarithmic slopes of Q + and Q − is generally smaller because the scaling of Q + with P z0 is weaker than the α model predicts due primarily to a weaker than expected scaling of H s with P z0 .

Comparison with ZEUS Results
Prior to this work, radiation MHD simulations performed in the radiation dominated regime have only been run with the ZEUS FLD algorithm (e.g. Hirose et al. 2009). In all cases, the simulations were reported to be consistent with thermal stability, in tension with our Athena results, for which thermal runaways are the norm. Hirose et al. (2009) explained this stability by noting that fluctuations in the magnetic field lead those in the radiation. And therefore they argue that magnetic field controls the radiation, not vice versa.
The picture outlined above also applies to the Athena simulations. In fact, we see an almost identical time lag between the fluctuations of radiation pressure and magnetic energy, with magnetic pressure leading the radiation pressure by ∼ 10 orbits. Therefore, differences in the stability properties of the two codes most likely originates in differences in the way the radiative cooling and the saturation level of the turbulent stresses respond to fluctuations in the pressure. The latter topic is still not well understood, even in isothermal stratified simulations, and requires a direct and detailed comparison between the codes. In collaboration with Hirose et al., we have initiated a detailed comparison between the codes and the simulations results, which will be the subject of a future paper (Jiang et al. in prep.). Since the focus of this paper is the Athena results, we simply summarize the most salient area of conflict or agreement between the simulations.
Overall, our solutions confirm many important results regarding the vertical structures and time evolution reported with ZEUS. As discussed in the previous sections, our simulations can run for many thermal times without any clear evidence of thermal runaway. In particular, we note that during the first 50 orbits of evolution in the moderately radiation dominated domains (RHVET and RMLVET), we find a relatively steady vertical structure is maintained for several (∼ 5) thermal times, characterized by large amplitude fluctuations in the stress and energy densities but with no secular trends. When a thermal runaway does eventually occur, it generally takes many thermal times for the expansion or contraction to proceed. Hence there is no sign of an exponentially growing fluctuations as one would expect from a linear instability.
The discrepancies between the results only becomes apparent on long time scales (i.e. many thermal times), since the Hirose et al. (2009) simulations do not show clear evidence of runaway behavior, even when evolved for much longer times (up to 600 orbits). The results of section 4.4 suggest that these discrepancies may be attributed, at least in part, to both the horizontal domain size and the use of different radiation transfer algorithms. It is reassuring that the simulations seem to agree when we use Athena's FLD algorithm with the same domain size and resolution as ZEUS.
In contrast, Athena simulations always end in a thermal runaway at the larger radiation-to-gas pressure ratios (∼ 200 in RSVET and RSFLD), independent of the transfer algorithm employed. Recently, experiments using ZEUS in the large radiation pressure, low surface density regime as in RSVET (Hirose et al., private communication 2013), find similar behavior. They perturbed the disk temperature by 2% to 10% after restarting evolution of a ZEUS simulation with the same parameters as RSVET at 100 orbits. Depending on the type of the perturbations they added, the disk expanded, collapsed, or continued running in an apparently stable equilibrium. When they continued to run the simulation to 300 orbits without adding any perturbation, the disk still collapsed. Therefore, it seems possible for ZEUS to find similar thermal runaways to those described above when the domains are very radiation dominated.

CONCLUSIONS
The evolution of MRI-unstable radiation pressure dominated accretion disks computed with the new radiationtransfer module in Athena always show thermal runaway. The disk can either expand or collapse, depending on the surface density of the disk, and numerical parameters of the calculations. The dependence of the heating and cooling rates, Q + and Q − respectively, on mid-plane pressure is consistent with the general criterion for thermal instability, in that Q + increases with mid-plane pressure faster than Q − .
For the strongly radiation pressure dominated case with lower surface density, the disk always undergoes a thermal runaway, independent of the radiation transfer algorithm and the numerical parameters. The time lag between pressure and magnetic energy density fluctuations, which is used to explain the results in Hirose et al. (2009), still exists in all these unstable solutions. This suggests that the time lag itself is not sufficient to stabilize the disk, which is consistent with the linear analysis including the time lag (Ciesielski et al. 2012).
For the marginally radiation pressure dominated case, the ZEUS code does not show thermal runaway behavior reported here even over much longer integration times (up to 600 orbits). To investigate this difference, we have implemented FLD in Athena. We are able to reproduce the behavior reported by Hirose et al. (2009), only if we use FLD and the small horizontal box size (∼ L x = 0.5H, L y = 2H) in Athena. When we use different radiation transfer algorithms for this small box size, or we increase the horizontal box size in Athena, we always find a thermal runaway. Therefore, the previous reports of thermal stability in the moderately radiation pressure dominated disk may be because of, in part, the smaller box size used. We believe simulations run on larger domains to be more robust (e.g., Simon et al. 2012).
A natural question about the observed expansion and collapse in our simulations is whether this behavior will continue if a larger vertical box size, or higher numerical resolution at the mid-plane, are used. For example, in those cases where the disk collapses, we cannot resolve the MRI modes that may still fit into the disk, and it is possible these modes might eventually lead to turbulence, heating, and re-expansion of the disk. However, in many respects this issue is moot. The fact that the mid-plane pressure can change by a factor of ∼ 10 has important consequences for the structure, evolution, and observational appearance of radiation dominated disks.
Moreover, if the vertical thickness is varying by such large amplitudes in local patches of the disk, we expect radial fluxes of mass, momentum, and energy may become important, which invalidates the use of the shearing box approximation. Thus, we consider it appropriate to characterize fluctuations of such large amplitude as a runaway, and moreover the final state reached in such runaways is clearly beyond the scope of local shearing box simulations, and requires global models.
All the simulations reported here are done with net azimuthal magnetic field. It would be very interesting to see how the radiation dominated disks behave with net vertical magnetic flux. MRI with net vertical magnetic flux can generate more vigorous turbulence than the net azimuthal magnetic field case and the stress increases with vertical box size in the unstratified simulations (Hawley et al. 1995). Therefore, we expect the stress will also increase with disk scale height and thermal instability should still exist in this case. However, future simulations are required to confirm this supposition The thermal runaways we find here may have important implications for observations. The α disk model predicts that the classical thermal instability will lead to a limit-cycle behavior, where the disk switches between a low temperature, gas pressure dominated state and high temperature, radiation pressure dominated state (Janiuk et al. 2002;Done et al. 2007). However, observations of most Galactic X-ray binaries do not find variability similar to that predicted (Gierliński and Done 2004), except for the well known source GRS1915+105 (Done et al. 2004a(Done et al. , 2007 and recently reported source IGR J17091-3624 (Altamirano et al. 2011). It is believed that the accretion rate reaches the Eddington limit in GRS1915+105, and this is the trigger for instability in this case (Done et al. 2004a,b). Our new simulations show that although thermal runaway may occur, it does so on growth rates that are smaller by a factor of a few compared to the α model. Moreover, when the radiation pressure is only a few times the gas pressure (as in RHVET and RMLVET), the disk can survive for several thermal times before runaway occurs. In this case the classical limit-cycle behavior will likely be modified. On the other hand, when the radiation pressure is hundred times the gas pressure (as in RSVET), the disk collapses within ∼ 2 thermal times. In this case, we expect that the classical limit-cycle behavior may emerge (Janiuk et al. 2002), which may have relevance to the observed variability of GRS1915+105.
Finally, we remark that the observational appearance of real disks will depend on how the thermal runaway saturates (Dexter and Quataert 2012). This cannot be addressed in the local shearing box approximation. Moreover, the thermal instability can be modified by radial advection of energy (Abramowicz et al. 1988), and by fluctuations in the surface density on long wavelengths (Lightman and Eardley 1974). All of these issues require global simulations of the saturation of the MRI in radiation dominated disks, which is the focus of our current effort.