RADIATION MAGNETOHYDRODYNAMIC SIMULATIONS OF THE FORMATION OF HOT ACCRETION DISK CORONAE

The standard α-disk model (Shakura & Sunyaev 1973) is usually used to account for the thermal emission from both active galactic nucleus (AGN; e.g., Krolik 1999) and X-ray binaries (e.g., Done et al. 2007). However, a hard X-ray tail, or a soft X-ray component with energies much larger than expected from thermal emission from the disk are often observed from both AGNs (Elvis et al. 1978) and X-ray binaries (Remillard & McClintock 2006; Done et al. 2007). During a state transition in X-ray binaries, the hard X-ray tail may actually become stronger than the thermal emission from the accretion disk (Remillard & McClintock 2006).

Hard X-rays from black hole accretion flows are usually believed to be formed via inverse Compton scattering of seed photons from the accretion disk by a hot corona above the surface (Bisnovatyi-Kogan & Blinnikov 1976; Haardt & Maraschi 1991, 1993; Svensson & Zdziarski 1994; Zdziarski et al. 1999). The corona, once formed, is usually thought to be located in a compact region near the innermost stable circular orbit (e.g., Reis & Miller 2013; Dexter & Blaes 2014). However, despite the rich observational evidence for the existence of hot coronae, there is no general consensus regarding the mechanism by which such coronae are formed. Many models (Haardt & Maraschi 1991; Svensson & Zdziarski 1994) simply assume that some fraction of the energy liberated by accretion is dissipated in optically thin regions. Magnetic reconnection is often invoked as the energy dissipation mechanism in the corona (e.g., Galeev et al. 1979; Goodman & Uzdensky 2008; Uzdensky & Goodman 2008). On the other hand, magneto-rotational instability (MRI) is now understood to be the mechanism that drives angular momentum transport and energy dissipation in black hole accretion disks (Balbus & Hawley 1991, 1998). It is therefore of great interest to explore whether the MRI naturally leads to high-temperature, magnetic pressure supported corona, and if so, by what mechanism.

In fact, early numerical simulations that studied the vertical structure of accretion disks with MRI turbulence (Miller & Stone 2000) showed that indeed a strongly magnetized region was formed in the low density upper layers of the disk. However, these simulations adopted an isothermal equation of state, and therefore, the gas temperature in the upper layers could only be estimated ex post facto. Using the flux-limited diffusion module (Turner & Stone 2001) in the ZEUS code, Hirose et al. (2006) studied the vertical temperature structure of accretion disks with MRI turbulence and radiative cooling. Although they also found magnetic pressure supported low density regions emerged in the upper layers of the disk, the gas temperature in these layers was always much smaller than at the disk mid-plane. Detailed spectral modeling (Blaes et al. 2006) based on the vertical profiles of dissipation and magnetic support taken from these simulations showed a slight hardening of the spectrum compared to previous models that neglected these effects, but no sign of a much hotter corona. Subsequent models computed using the same code (Krolik et al. 2007; Blaes et al. 2007, 2011; Hirose et al. 2009b) that explored a wide range of ratios between radiation and gas pressure always found similar temperature profiles. At the same time, a common property of these simulations is a significant increase in dissipation rate per unit mass in the surface regions, suggesting a temperature inversion might be possible under the right circumstances.

All previously published radiation magnetohydrodynamic (MHD) simulations of the MRI have considered a large surface density so that the total electron scattering optical depth from the disk mid-plane to the surface is ∼10⁴. In this paper, we report on new radiation MHD simulations performed using the full-transport radiation transfer module in Athena Davis et al. (2012) and Jiang et al. (2012). We find that when the disk surface density is decreased so that the total electron scattering optical depth is only ∼10², a strong temperature inversion forms above the disk surface due to dissipation of the turbulence in magnetically supported regions which are optically thin. This temperature inversion is (in many respects) consistent with the observationally inferred properties of accretion disk coronae, and we will refer to it as a corona throughout this work. Because of the thermal runaways that we generically observe in radiation pressure dominated simulations (Jiang et al. 2013a), we focus solely on gas pressure dominated disks in this paper. This constraint prevents the simulations from producing the very high temperatures (∼10⁹ K) that are observationally inferred but does not prevent an initial exploration of this potentially new mechanism for corona formation.

We adopt the local shearing box approximation (Hawley et al. 1995; Jiang et al. 2013a), which means we study a local patch of the accretion disk in a frame rotating with orbital frequency Ω at a fiducial radius r₀ from the central black hole (BH) with mass M_BH = 6.62 M_☉. Curvature of the orbit is neglected so the radial, azimuthal and vertical directions are represented by the local Cartesian coordinate (x, y, z) with unit vector ($\hat {\boldsymbol {i}}$, $\hat {\boldsymbol {j}}$, $\hat {\boldsymbol {k}}$), respectively. The vertical component of the gravitational force from the BH under the thin disk approximation is included as $-\rho \Omega ^2z \hat{\boldsymbol {k}}$. Energy changes of the gas due to Compton scattering is approximated as (Hirose et al. 2009b) ΔE_s = −4cE_rρκ_es(T − T_r)/T_e, where c is the speed of light, E_r is the radiation energy density, ρ is gas density, and electron scattering opacity κ_es = 0.33 cm² g⁻¹, T is the gas temperature, and T_r is the radiation temperature defined as T_r ≡ (E_r/a_r)^1/4 with radiation constant a_r = 7.57 × 10¹⁵ erg cm⁻³ K⁻⁴. The equivalent electron temperature is defined as T_e ≡ m_ec²/k_B = 5.94 × 10⁹ K, where m_e is the electron mass and k_B is the Boltzmann constant. The radiation field is always assumed to be a Planck distribution in this formula for Compton scattering.

The complete set of radiation MHD evolutionary equations we solve are given by Equations (2) and (3) of Jiang et al. (2013a). Our solutions are computed using the Godunov radiation MHD code based on a variable Eddington tensor (VET) as described and tested by Jiang et al. (2012) and Davis et al. (2012), with the improvements described in the Appendix of Jiang et al. (2013b). The adiabatic index is chosen to be γ = 5/3 with mean molecular weight 0.6. Plank-mean free–free absorption opacity $\kappa _{aP}=3.7\times 10^{53}(\rho ^9/E_g^7)^{1/2}$ cm² g⁻¹ and Rosseland-mean free–free absorption opacity $\kappa _{aF}=1.0\times 10^{52}(\rho ^9/E_g^7)^{1/2}$ cm² g⁻¹, where E_g is the gas internal energy density. The unit of the magnetic field is chosen such that magnetic permeability is one (Stone et al. 2008).

2.1. Initial and Boundary Conditions

During the linear growth phase of the MRI in the first ∼10 orbits of both simulations, the disks cool radiatively through the surfaces and slowly collapse. After the MRI saturates, the disks are heated by the dissipation of turbulence. Subsequently, the disks adjust their vertical structure to reach a thermal equilibrium state in which radiative cooling and turbulent heating are balanced. Since the disks are supported by gas pressure, they are thermally stable (Shakura & Sunyaev 1976; Jiang et al. 2013a). We run simulation A for more than 120 orbits with average cooling time to be five orbits. Simulation B runs to 200 orbits, much longer than the average cooling time of 10 orbits. Here the cooling time is defined to be the ratio between the total energy density inside the simulation domain and the cooling rate from the surfaces of the domain. If we map the local shearing box simulations to global disk structures by assuming cooling is balanced by the release of gravitational energy, the steady state of simulation A roughly corresponds to a thin disk model with α = 0.03 and 3.8 × 10⁻⁵ Eddington accretion rate, while simulation B corresponds to a thin disk model with the same α but 0.12 Eddington accretion rate. That is why although simulation A is closer to the central black hole, it can have a smaller surface density and lower mid-plane temperature.

The evolution histories for density ρ, temperature T, and azimuthal magnetic field B_y for simulation A are shown in a space–time diagram in Figure 1, generated by averaging quantities across horizontal planes and then plotting the resulting vertical profiles versus time. Consistent with previous simulations with much larger surface density (Hirose et al. 2006; Jiang et al. 2013a), the sign of B_y flips roughly every 10 orbits, which is the well-known butterfly diagram observed in almost all the vertically stratified shearing box simulations (e.g., Stone et al. 1996; Miller & Stone 2000; Davis et al. 2010; Jiang et al. 2013a). The space–time diagram for density also shows very similar structures. As the magnetic field rises buoyantly, it carries along some gas which is denser than the surrounding coronal material; this gas later falls back toward the disk mid-plane. During the whole simulation, the photosphere is well inside the simulation box, as labeled by the white line at the top panel of Figure 1. Simulation B also shows very similar behavior.

Zoom In Zoom Out Reset image size

The most dramatic difference between simulations A and B is the space–time diagram of the gas temperature. Simulation B is similar to previous MRI simulations with radiation transfer (Hirose et al. 2006) in that gas temperature peaks at the disk mid-plane and drops with height. However, for simulation A, as shown in Figure 1, the gas temperature above the electron scattering photosphere is dramatically larger than the gas temperature at the disk mid-plane. The gas temperature is maximum during the periods when the magnetic field rises buoyantly. This correlation exists because dissipation in the corona is directly from the magnetic energy density, which is not amplified locally at each height but carried here from the disk mid-plane.

3.1. Vertical Profiles of the Disk

To investigate the vertical structure of the disk quantitatively, we time-average the vertical profiles between 25 and 120 orbits for simulation A and 50–200 orbits for simulation B. The resulting time-averaged vertical profiles of gas (T) and radiation (T_r) temperature, gas (P_g), radiation (P_r) and magnetic pressure (P_m), as well as three components of the Eddington tensor (f_xx, f_yy, f_zz) are shown in Figure 2. Within the disk photosphere, gas and radiation temperature are coupled, and they both decrease with height from the disk mid-plane. In this optically thick region, the Eddington tensor is close to $1/3{{\sf \boldsymbol I}}$ as expected and P_g is much larger than P_m and P_r. Once above the photosphere, T and T_r decouple. For simulation B, T continues to decrease with height. The rapid change near the boundary is likely caused by the boundary condition. The whole temperature profile is consistent with Figure 3 of Hirose et al. (2006). However, for simulation A, radiation temperature T_r stays almost flat above the photosphere while gas temperature increases very quickly by a factor of ∼30. Gas pressure also drops quicker than magnetic pressure above the photosphere, and this region is supported primarily by magnetic pressure. High gas temperature with strong magnetic pressure support are the defining characteristics of the corona observed in our simulations. The sharp increase of the gas temperature also causes peaks in the horizontal components of the Eddington tensor and a dip in the vertical component. This is because the rising temperature in combination with the falling density produces a localized (in z) maximum in the emissivity. At nearly horizontal viewing angles, one has a longer line of sight through this emissivity "bump" yielding larger intensities for horizontal rays and lower intensities for vertical rays. This projection effect becomes less important closer to the top of the domain and the intensity becomes more nearly isotropic and eventually slightly limb-darkened. After the gas temperature becomes flat, f_zz becomes larger than 1/3 while f_xx and f_yy becomes smaller than 1/3, which are consistent with profiles of Eddington tensor in simulation B.

Zoom In Zoom Out Reset image size

3.2. The Dissipation Profile

The reason that gas temperature increases above the photosphere is that there is still significant dissipation in this region. Because it is optically thin, the gas cannot cool easily, and its temperature rises. Because of the different temperature profiles between simulation A and B, we naturally expect different dissipation profiles for the two simulations with significantly different surface densities. For local shearing box simulations of MRI, the energy input is proportional to the sum of Maxwell and Reynolds stress W_xy from MRI turbulence (Hawley et al. 1995; Jiang et al. 2013b). The vertical profiles of dissipation can therefore be directly measured from our simulations, removing one of the largest uncertainties inherent in emission models of accretion disks (Davis et al. 2005; Tao & Blaes 2013). The horizontally and time averaged profiles of stress per unit mass as a function of total optical depth τ and surface density σ for simulation A and B are shown in Figure 3. As this figure shows, W_xy/ρ increases rapidly with height, while the standard α disk model assumes it is a constant. By comparing simulation A and B, it is clear that the vertical profile of W_xy/ρ becomes steeper as surface density decreases. At the same time, the position of the photosphere moves closer to the disk mid-plane. In simulation B, only 0.0094% of the total mass is located above the photosphere, which only contains 0.085% of the energy dissipation. However, for simulation A with much smaller surface density, the optically thin photosphere includes 0.26% of the total mass and 3.4% of the total dissipation. Clearly, a larger fraction of energy is dissipated in the optically thin region in simulation A, and ultimately, this is why a hot corona is formed. This also suggests that the typical scale height of dissipation does not scale with density scale height. This is not surprising as dissipation in MRI is ultimately related to the magnetic field, which always tends to rise up to the low density region due to buoyancy.

Zoom In Zoom Out Reset image size

To provide some context for how dissipation occurs in the simulation, Figure 4 shows a snapshot of the typical magnetic field structure in the top half of the box. The magnetic field is mainly along the toriodal direction in the optically thick mid-plane, where the Maxwell stress is the largest. In the magnetic pressure supported optically thin photosphere, a significant vertical component of the magnetic field can be formed due to Parker instability (Blaes et al. 2007; Hirose & Turner 2011). To investigate the possible association of these structure with dissipation, we plot the time averaged vertical profiles of the current density J, J², and their standard deviations normalized to horizontally averaged values in Figure 5. We expect that regions where the fluctuations in these quantities are large are likely to be dominated by dissipation at small scales. The figure demonstrates the horizontal fluctuations are larger in the corona region compared with the disk mid-plane for both J and J², which suggests that small-scale dissipation is still dominant in the corona. We caution, however, that it is not clear if fluctuations in J or J² are indeed a good indicator of the scales in which dissipation occurs, and moreover, this dissipation may not necessarily be related to reconnection. A more quantitative exploration of the role of reconnection in MRI turbulence in both the mid-plane and corona is beyond the scope of this work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By comparing two different simulations with MRI turbulence and radiation transfer, we first confirm that when the surface density is so large that the total electron scattering optical depth is ∼10⁴, there is very little dissipation in the optically thin part, and the temperature in the magnetic pressure supported region is always much smaller than the disk mid-plane. When the surface density of the disk is small enough such that the total electron scattering optical depth is only ∼200, a non-negligible fraction of dissipation from MRI occurs above the photosphere, and a high temperature, magnetic pressure supported corona is self-consistently generated. Note that the disk is still in the gas pressure dominated regime, and the temperature in the corona is only ∼10⁷ K in this particular run. Since this is still much smaller than what is required to make the observed hard X-rays, we do not claim to have reproduced the corona of observed systems. Observationally consistent coronal models for most observed systems will require a larger fraction of the dissipation above the photosphere than found in simulation A (∼3.4%). Nevertheless, we consider this a useful demonstration that MRI turbulence can produce temperature inversions in magnetically dominated surface regions in simulations with realistic thermodynamics. This also motivates future studies to see how the structures of the accretion disks will be changed if more and more energy is dissipated in the optically thin part when the disk surface density decreases further.

Figure 2 shows a notable contrast between the volume average and density-weighted average of the temperature in the surface regions. Near the surface, there is a significant anti-correlation between the temperature and density that leads to the volume average temperature being higher than the density-weighted average. Dissipation in this region is dominated by work done by magnetic stresses and numerical reconnection. The associated heating is rather inhomogeneous as indicated by Figure 5 and tends to deposit a larger amount of energy per unit mass in lower density regions, giving rise to the temperature-density anti-correlation. Such "patchy" heating is interesting because it may be necessary to explain the hard X-ray tails observed in high/soft or very high state X-ray binaries (see e.g., Section 5.3.1 of Done et al. 2007). In these systems, the hard spectral slopes are suggestive of high temperatures (or high electron energies in a non-thermal electron distribution), but the fraction of the power emitted in the component suggests that only a small fraction of photospheric photons are scattered. Such a geometry arises naturally if a significant fraction of the heating occurs in a small fraction of the volume, as the simulations imply.

One limitation of these shearing box simulations is that they consider only a small patch of the disk, so the global structure of the accretion disk and corona cannot be studied. Postprocessing of global MHD simulations with Monte Carlo radiative transfer calculations suggests that many aspects of the observed hard X-ray spectrum can be reproduced (Schnittman et al. 2013). One drawback of these calculations is that they rely on a heuristic cooling function to maintain a thin disk in the simulations (Noble et al. 2009). In the future, we hope to perform a similar analysis with global radiation MHD simulations, particularly since many of the observed systems are in the radiation pressure dominated regime. Here, we simply discuss how these local patches of the disk would fit into standard accretion disk models.

In the α disk model, the surface density of the disk decreases as the disk's effective and mid-plane temperatures decrease in the gas pressure dominated regime (Shakura & Sunyaev 1973), which is consistent with both ZEUS and Athena simulation results (Hirose et al. 2009a; Jiang et al. 2013a). In a global accretion disk model, this corresponds to reducing the accretion rate and total luminosity. But at the same time, we have shown that a larger fraction of energy is dissipated in the corona region. Therefore, the ratio between luminosity from the high energy band and the low energy band, or the spectral index, would increase when the total luminosity decreases. The transition from the high–soft state to the low–hard state in X-ray binaries is a case where spectrum hardens in this manner as the luminosity decreases (Remillard & McClintock 2006). Even if the inner disk is radiation dominated at the beginning of the transition, it should become gas pressure dominated as the accretion rate drops so the surface density must decrease as well. There is some observational evidence that during this transition, the optically thick disk is truncated at radii larger than the last stable circular orbit (Done 2010 and references therein). In the inner disk, the geometrically thin and optically thick disk must become a geometrically thick and optically thin advection-dominated accretion flow (ADAF)-like flow (Narayan & Yi 1994). If this interpretation is correct, there is likely to be some transition region with a total electron scattering optical depth smaller than ∼100, which corresponds to the parameters adopted in our simulations.

For the radiation pressure dominated regime, the α disk model predicts that the disk effective and midplane temperatures will increase as the disk surface density decreases. However, our shearing box simulations of the radiation pressure dominated regime show thermal runaways (Jiang et al. 2013a), and it remains unclear how these runaways might saturate in a global disk. Nevertheless, Jiang et al. (2013a) found that in the radiation pressure dominated regime, the dissipation scale height changes slower than the change of density scale height. Therefore, a larger fraction of energy will be dissipated at high altitude when disk surface density decreases. On one hand, this may be inconsistent with some high/soft state observations which have relatively low coronal contribution. On the other hand, the fact that the very high state (or steep power law state) has a strong tail and only seems to occur at high luminosities (Remillard & McClintock 2006) suggests that the surface density of the accretion disk is only small enough to enable the formation of corona when the accretion rate is very large in the radiation pressure dominated regime. Therefore, it will still be interesting to see in the radiation pressure dominated case before the thermal runaway completely changes the structures of the disk, how corona will be formed when the surface density is decreased. This will be studied with future shearing box simulations.

Most AGNs are expected to be in the radiation pressure dominated regime, except when at low Eddington rates where a transition to an ADAF-like flow occurs. Therefore, we can only speculate on the observational implications of our current, gas pressure dominated results. We note that X-ray and optical observations of AGNs do indicate a trend of hardening of the spectral index as luminosity decreases (Steffen et al. 2006; Just et al. 2007), although the black hole mass is uncertain so it is unclear what the relevant Eddington ratios are in these systems. Recently, a new population of quasars was found to show much weaker X-ray emission compared with other AGNs with similar luminosity range (Wu et al. 2011; Luo et al. 2013). Observationally, it is still debated whether these quasars are intrinsically X-ray weak or obscured. If the AGN X-ray emission is indeed from the corona that is generated by MRI turbulence as found here, then these quasars may be intrinsically X-ray weak when most of the dissipation happens inside the photosphere and when the corona becomes very weak or disappears. Finally, we note that the increased dissipation per unit mass near the disk surface may be the mechanism needed to generate the continuum of soft X-ray excess inferred in some of the AGNs (Czerny et al. 2003; Done et al. 2012). Ultimately, testing these hypotheses will require global simulations of accretion flows in both the gas and radiation pressure dominated regimes.

Y.F.J. thanks Jenny Greene, Yue Shen, and Jianfeng Wu for helpful discussions on the X-ray and optical observations of AGNs, and J. Goodman for discussions which motivated this project. We also thank the anonymous referee for helpful comments that improved the paper. This work was supported by the NASA ATP program through grant NNX11AF49G, and by computational resources provided by the Princeton Institute for Computational Science and Engineering. Some of the simulations were performed on the Pleiades Supercomputer provided by NASA. Y.F.J. is supported by NASA through Einstein Postdoctoral Fellowship grant number PF-140109 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. This work was also supported in part by the U.S. National Science Foundation, grant NSF-OCI-108849 and NSF-AST-1333091.