Global Radiation Magneto-hydrodynamic Simulations of Sub-Eddington Accretion Disks around Supermassive Black Holes

TL;DR

This study employs 3D radiation magneto-hydrodynamic simulations to explore the structure and behavior of sub-Eddington accretion disks around supermassive black holes, revealing differences from standard models and implications for AGN coronae.

Contribution

It provides new insights into the magnetic and radiation-driven dynamics of accretion disks at sub-Eddington rates, challenging standard thin disk assumptions.

Findings

Disks are supported by magnetic pressure with larger scale height.

No thermal instability observed over many thermal timescales.

Coronae are generated only in the inner regions, with sizes consistent with observations.

Abstract

We use global three dimensional radiation magneto-hydrodynamic simulations to study the properties of inner regions of accretion disks around a 5\times 10^8 solar mass black hole with mass accretion rates reaching 7% and 20% of the Eddington value. This region of the disk is supported by magnetic pressure with surface density significantly smaller than the values predicted by the standard thin disk model but with a much larger disk scale height. The disks do not show any sign of thermal instability over many thermal time scales. More than half of the accretion is driven by radiation viscosity in the optically thin corona region for the lower accretion rate case, while accretion in the optically thick part of the disk is driven by the Maxwell and Reynolds stresses from MRI turbulence. Coronae with gas temperatures > 10^8 K are generated only in the inner \approx 10 gravitational radii in both simulations, being more compact in the higher accretion rate case. In contrast to the thin disk model, surface density increases with increasing mass accretion rate, which causes less dissipation in the optically thin region and a relatively weaker corona. The simulation results may explain the formation of X-ray coronae in Active Galactic Nuclei (AGNs), the compact size of such coronae, and the observed trend of optical to X-ray luminosity with Eddington ratio for many AGNs.

Tables

Table 1. Table 1 Simulation Parameters

Variables/Units	AGN0.2	AGN0.07
$r_i/r_g$	80	80
${\rho }_i/{\rho }_0$	1	1
$T_i/T_0$	3.18	3.18
$⟨P_r/P_g⟩$	$4.47\times {10}^5$	$4.61\times {10}^5$
${⟨P_r/P_g⟩}_{\rho }$	$4.40\times {10}^2$	$4.39\times {10}^2$
$⟨P_B/P_g⟩$	$1.27\times {10}^{-2}$	$2.32\times {10}^{-4}$
${⟨P_B/P_g⟩}_{\rho }$	$7.70\times {10}^{-3}$	$7.32\times {10}^{-5}$
$\Delta r/r$	$6.1\times {10}^{-3}$	$6.1\times {10}^{-3}$
$\Delta \theta$	$6.1\times {10}^{-3}$	$6.1\times {10}^{-3}$
$\Delta \varphi$	$6.1\times {10}^{-3}$	$6.1\times {10}^{-3}$
$N_n$	80	80

Equations

$$P_{r,\mathrm{vis}}^{r\varphi }=-\frac{8}{27}\frac{E_r}{\rho {\kappa }_{R}c}{D}_{r\varphi },$$

$$P_{r,\mathrm{vis}}^{r\varphi }=-\frac{8E_r}{27c}k(\tau )D_{r\varphi }^{\prime },$$

Full text

1. Introduction

Active Galactic Nuclei (AGNs) are believed to be powered by accretion onto the central black holes. Luminosity of most AGNs is observed to be smaller than the Eddington limit as defined by the electron scattering opacity $L_{\rm{Edd}}=1.5\times 10^{46}M_{\rm{BH}}/\left(10^{8}M_{\odot}\right)$ erg s*-1*. The standard thin disk model (Shakura & Sunyaev 1973) is usually adopted to describe the accretion disks in this regime when the AGN luminosity is larger than $\approx 0.01L_{\rm{Edd}}$. The key assumptions in this model are that the disk is optically thick and supported by thermal pressure (gas or radiation), and all the dissipation is radiated away locally. The detailed structures of the disk are then determined by the famous $\alpha$ assumption, which assumes that the stress responsible for the angular momentum transport in the disk is proportional to the total pressure with a constant value.

Comparing the predictions from the standard thin model with the observed properties of AGNs has raised many questions on the assumptions in this model (Koratkar & Blaes 1999). The standard thin disk model subjects to the well known thermal and inflow instabilities (Shakura & Sunyaev 1976; Lightman & Eardley 1974; Piran 1978), which should cause the disk to evolve away from the equilibrium state in a few thermal time scales (Jiang et al. 2013a; Fragile et al. 2018) and probably show some kind of limit-cycle behavior (Honma et al. 1991; Janiuk et al. 2002). However, this is not observed for most AGNs. The predicted spectrum from this model is also not consistent with the observed energy distribution (Zheng et al. 1997; Davis et al. 2007; Laor & Davis 2014). Particularly, the predicted edge feature in the spectrum is also not observed (Sincell & Krolik 1997; Shull et al. 2012; Tilton et al. 2016).

The standard thin disk model assumes that thermal pressure (gas or radiation) is dominant in the disk. Alternative possibility that the vertical component of gravity in the accretion disk may be primarily balanced by magnetic pressure has been proposed (Shibata et al. 1990; Pariev et al. 2003). Begelman & Pringle (2007) studied the structures and stability of magnetic pressure dominated disks by assuming the toroidal magnetic fields will saturate to a level so that the associated Alfvén velocity reaches $\sqrt{c_{s}V_{k}}$ (Pessah & Psaltis 2005; Das et al. 2018), where $c_{s}$ is the gas sound speed and $V_{k}$ is the disk rotation speed. Amplification of magnetic fields near the disk midplane is thought to be balanced by the escape of magnetic fields away from the midplane due to buoyancy. This magnetic elevated model is found to have a larger pressure scale height and does not subject to the thermal and viscous instabilities (Sa̧dowski 2016), which have interesting implications for both X-ray binaries and AGNs (Begelman et al. 2015; Begelman & Silk 2017; Dexter & Begelman 2019). However, assumptions in the magnetic elevated disk model have not been checked numerically. It is also unclear how the strong radiation pressure in AGN accretion disks will modify the structures, which is typically neglected in these models. It is known that in order to reach the magnetic pressure dominated state, large scale poloidal or radial magnetic fields are required (Salvesen et al. 2016b, a; Fragile & Sa̧dowski 2017). The advantage of our simulations is that gas, radiation and magnetic pressure are all included self-consistently. Our goal is to study the structures of the disk that will be formed with these magnetic field configurations for realistic parameters of AGNs. Our simulations form the a disk in a region that is initially vacuum by accreting gas from larger radius. Although the resulting accretion disk depends on our initial condition for the magnetic fields, the disk forms self-consistently. This distinguishes it from previous work that initialized the simulation with a standard thin disk model and study its long term evolution (Fragile et al. 2018).

Formation of coronae in AGNs, as well as the dependence of corona properties on the accretion flow, remain a puzzle. The effective temperature of the accretion disks in AGNs (typically $\sim 10^{4}-10^{5}$ K) is too low to produce the widely observed X-rays in AGNs via thermal emission. It is proposed that high temperature corona, as inspired by the solar corona, is formed on top of the accretion disk (Bisnovatyi-Kogan & Blinnikov 1976; Haardt & Maraschi 1991, 1993; Svensson & Zdziarski 1994; Zdziarski et al. 1999). X-rays are produced via Compton scattering of the seed photons emitted by the disk with the hot electrons in the corona, which is infereed to exist in a compact region near the black hole ($\sim 10$ gravitational radii, Reis & Miller 2013; Uttley et al. 2014). The amount of X-rays that can be produced via this mechanism compared with the thermal emission from the disk is determined by the fraction of accretion power that is dissipated in the corona region. This is basically a free parameter in these models since the mechanism to produce the corona in AGNs is unknown. Since the early isothermal simulations of vertically stratified accretion disk with turbulence generated by the Magneto-rotational instability (MRI) (Miller & Stone 2000), it is commonly observed that magnetic fields amplified by MRI near the midplane of the disk can buoyantly rise to the surface and create a magnetic pressure dominated low density region, which is suspected to be the corona. However, most local simulations that determine the thermal properties of the disks based on radiative cooling (using both the flux-limited diffusion and VET method) typically find that the gas in the magnetic pressure dominated region is not heated to a very high temperature due to insignificant dissipation (Krolik et al. 2007; Blaes et al. 2007, 2011; Hirose et al. 2009; Jiang et al. 2013a, 2016). An exception is Jiang et al. (2014b), which shows that the amount of dissipation in the magnetic pressure dominated region can be increased by reducing the surface density of the disk, which can increase the gas temperature in this region. Surface density is a free parameter in these local shearing box simulations and the value required to produce the high temperature corona is smaller than what the thin disk model will predict for the same accretion rate. It is therefore necessary to use global simulations to determine the surface density of the disk as well as the properties of coronae self-consistently for a given mass accretion rate, which is one goal of the paper.

Jiang et al. (2016) shows that for typical density and temperature expected for AGN accretion disks based on the standard thin disk model, the Rosseland mean opacity should be larger than the electron scattering value due to irons. The density and temperature dependences of the iron opacity peak can modify the thermal stability and structures of AGN accretion disks significantly. However, the local shearing box simulations done by Jiang et al. (2016) adopt the surface density as given by the standard thin disk model. Whether this opacity peak will show up or not depends on the actual disk structure we will get. We will include the full opacity table in the simulations to capture its potential importance.

The remainder of this paper is organized as follows. In Section 2, we describe the simulation setup. Detailed structures of the disk are described in Section 3. We discuss the implications of our simulations in Section 4.

2. Simulation Setup

We solve the same set of ideal MHD equations coupled with the time dependent radiative transfer equation for specific intensities as in Jiang et al. (2017a) using the code Athena++ (Stone et al, in preparation). We carry out two simulations, AGN0.2 and AGN0.07, for a $M_{\rm{BH}}=5\times 10^{8}M_{\odot}$ black hole with accretion rates smaller than the Eddington value $\dot{M}_{\rm{Edd}}\equiv 10L_{\rm{Edd}}/c^{2}=8.22\times 10^{26}\ {\rm g/s}$. The simulations are performed with the pseudo-Newtonian potential (Paczyńsky & Wiita, 1980) $\phi=-GM_{\rm{BH}}/(r-2r_{g})$ to mimic the general relativity effects around a Schwarzschild black hole, where $G$ is the gravitational constant while gravitational radius $r_{g}\equiv GM_{\rm{BH}}/c^{2}=7.42\times 10^{13}\ {\rm cm}$. We initialize a torus centered at $80r_{g}$ with the maximum density $\rho_{0}=10^{-8}\ {\rm g}\ {\rm cm}^{-3}$ and temperature $3.18T_{0}$, where $T_{0}=2\times 10^{5}{\rm K}$. The shape of the torus is the same as the ones used in Jiang et al. (2017a). The inner edge of the torus is at $40r_{g}$ and the region inside that radius is filled with density floor $10^{-9}\rho_{0}$ initially. The initial properties of the torus, including the ratios between the averaged radiation pressure $P_{r}$, gas pressure $P_{g}$ and magnetic pressure $P_{B}$, are summarized in Table 1. The main difference between the two simulations is the initial magnetic field in the torus. The run AGN0.2 uses a single loop of poloidal magnetic field while the run AGN0.07 adopts multiple magnetic field loops as described in Jiang et al. (2017a). The different setups lead to different mass accretion rates in the disks that are formed near the black hole.

We use four levels of static mesh refinement to cover the whole simulation domain $\left(r,\theta,\phi\right)\in\left(4r_{g},1600r_{g}\right)\times\left(0,\pi\right)\times\left(0,2\pi\right)$. The level with the highest resolution reaches $\Delta r/r=\Delta\theta=\Delta\phi=6.1\times 10^{-3}$ for the region $\left(6r_{g},200r_{g}\right)\times\left(1.48,1.66\right)\times\left(0,2\pi\right)$, which covers most of the mass near the disk midplane. The equivalent resolution is $1024\times 512\times 1024$, which is necessary to resolve these sub-Eddington accretion disks.

We use $80$ discrete angles in each cell to resolve the angular distribution of the radiation field, which is in thermal equilibrium initially in the torus. We calculate the Rosseland mean opacity in each cell by using local density and temperature based on the OPAL opacity table with solar metallicity (Paxton et al. 2013; Jiang et al. 2015). Planck mean free-free absorption opacity is also included as in Jiang et al. (2017a). Each simulation takes $\approx 20-30$ millions CPU time in the ALCF machine Mira.

3. Results

3.1. Resolution for MRI

To quantify how well the MRI is resolved in our simulations, we calculate the quality factors $Q_{\theta}$ and $Q_{\phi}$, which are the ratios between the wavelength of the fastest growing MRI mode $\lambda=2\pi\sqrt{16/15}|v_{A}|/\Omega$ and cell sizes $r\Delta\theta,r\sin\theta\Delta\phi$ along $\theta$ and $\phi$ directions, where the Alfvén velocity $v_{A}$ is calculated for $B_{\theta}$ and $B_{\phi}$ respectively. Resolution studies for non-radiative ideal MHD simulations (Hawley et al., 2011; Sorathia et al., 2012) find that when $Q_{\phi}\gtrapprox 25,Q_{\theta}\gtrapprox 6$ or both $Q_{\phi}$ and $Q_{\theta}$ are larger than $10$, properties of MRI turbulence is converged with respect to resolution. Although there is no general criterion on the convergence of radiation MHD simulations of accretion disks, we use it as a way to compare our resolutions with these calculations. We calculate the azimuthally averaged quality factors during the final turbulent states for both of the calculations. For the simulation AGN0.07 at the disk midplane, $Q_{\phi}=100$ at $r=5r_{g}$ and decreases to 80 at $r=30r_{g}$, while $Q_{\theta}$ varies from 10 to 8 in the same radial range. For a fixed radius, the combined effects of increasing Alfvén velocity and reduced resolution with height from the disk midplane cause the quality factors to decrease by a factor of $\approx 2$ at the photosphere. For the simulation AGN0.2, the quality factors are larger as the disk is thicker. The midplane $Q_{\phi}$ changes from 200 at $5r_{g}$ to 40 at $40r_{g}$ while $Q_{\theta}$ changes from $\approx 20$ to $\approx 10$ between $5$ and $40r_{g}$ at the disk midplane. For this run, the quality factors increase by a factor of $\approx 2$ from the midplane to the photosphere. Therefore, the resolutions are sufficient to resolve the MRI turbulence according to the criterion found by non-radiative MHD simulations. This is possible for these sub-Eddington accretion disks because magnetic pressure is the dominant pressure, which makes the wavelength of the fastest growing MRI mode be comparable to the disk scale height and easier to resolve.

3.2. Simulation Histories

Gas in the initial torus flows towards the black hole and the accretion disk is slowly built up. Histories of the mass accretion rate $\dot{M}\equiv\int_{0}^{2\pi}\int_{0}^{\pi}\rho v_{r}r^{2}\sin\theta d\theta d\phi$ at $10r_{g}$ for the two runs are shown in the left panel of Figure 1. After a period of initial transient, the simulation AGN0.07 reaches an accretion rate $7\%\dot{M}_{\rm{Edd}}$ when averaged after time $t=2.4\times 10^{4}r_{g}/c$, which roughly corresponds to $13$ thermal time scales at $10r_{g}$. For the run AGN0.2, the averaged accretion rate is $22\%\dot{M}_{\rm{Edd}}$ after $t=3\times 10^{4}r_{g}/c$, which also corresponds to more than $\approx 6$ thermal time scales at $10r_{g}$. All the analysis for the time averaged properties of the two simulations are also performed for the same time intervals. The mass accretion rates show significant fluctuations due to turbulence. The standard deviation of $\dot{M}$ during the steady state for the run AGN0.07 is only $1.2\%$ of the averaged value while the corresponding value is $12.3\%$ for the run AGN0.2. The main difference of the two simulations is the magnetic field topology as described in Section 2. The run AGN0.2 adopts a single loop of magnetic field in the initial torus and there are net poloidal magnetic fields $\overline{B}_{\theta}$ through the inner region of the disk. The ratio between radiation pressure and $\overline{B}_{\theta}^{2}/2$ near the disk midplane is $\approx 10^{3}$. For the run AGN0.07, quadrupole magnetic fields are used in the initial condition, which result in net radial magnetic fields $\overline{B}_{r}$ near the disk midplane, although the shell averaged $B_{r}$ across the whole disk is still zero. The net $\overline{B}_{r}$ near the disk midplane is sheared into toroidal magnetic field within an orbital time scale and quickly builds up strong magnetic pressure, which then escapes from the midplane due to buoyancy. MRI turbulence is still developed in this case (Pessah & Psaltis 2005; Das et al. 2018), but it shows less variability compared with the other run. Detailed investigations of magnetic pressure dominated disks resulting from this magnetic field configuration are described in Section 3.3.

In most MHD simulations of accretion disks where radiative transfer is not calculated self-consistently, density is usually used as the proxy to understand the observed properties of these systems. Since our simulations calculate the photons emitted by the disk directly, we can calculate the frequency-integrated lightcurves for the time scale that our simulations can cover. In order to avoid contamination from photons generated by the torus at large radii where the disk has not reached steady state, we convert the radiation flux {{\mbox{\boldmathF}}_{r}} to the radial $F_{R}$ and vertical components $F_{z}$ in cylindrical coordinates and then integrate the total radiative luminosity leaving the cylindrical surface at $R=10r_{g}$ and height $z=100r_{g}$, which is well beyond the photosphere. The resulting lightcurves are shown in the right column of Figure 1, which have much weaker short time scale variabilities compared with what the variability in $\dot{M}$ might suggest. This is likely due to the scattering of photons through the optically thick disks. In particular, the power spectrum of the lightcurve from the lower accretion rate run (AGN0.07) has more power at high frequencies compared with the lightcurve power spectrum from the run with higher accretion rate (AGN0.2). This is because the midplane optical depth increases with increasing mass accretion rate (see Figure 6). The changes of variability amplitude with luminosity are also different in the two lightcurves, revealing different dynamo actions in the disk, which will be discussed in Section 4.

Histories of poloidal profiles of azimuthally averaged density $\rho$, gas temperature $T_{g}$ and toroidal magnetic field $B_{\phi}$ at $10r_{g}$ for the two runs are shown in Figure 2. The disk photosphere only extends to $\approx 6^{\circ}$ away from the disk midplane for the run AGN0.07 while in the simulation AGN0.2, the disk photosphere covers more than $14^{\circ}$. Once above the photosphere as indicated by the blue lines in the top panels, gas temperature increases rapidly to $\approx 10^{9}$ K, which is the corona region. More detailed properties of the corona are discussed in Section 3.6. In the optically thick part of the disk, gas and radiation are in thermal equilibrium with a temperature around $10^{5}$ K. The grid scale variation of gas temperature is unrealistic because gas pressure is smaller than $0.1\%$ of radiation and magnetic pressure and it suffers from numerical noise. The toroidal magnetic field switches signs near the disk midplane after $\approx 10^{4}r_{g}/c$ for the run AGN0.2, which is the well known butterfly diagram caused by the MRI dynamo (Stone et al., 1996; Miller & Stone, 2000; Davis et al., 2010; O’Neill et al., 2011; Simon et al., 2012; Jiang et al., 2013a, 2014a). This does not happen for the run AGN0.07, although MRI turbulence has also developed there. This is because shearing of the net radial magnetic field near the disk midplane for the run AGN0.07 always forms $B_{\phi}$ with the same sign in addition to the $B_{\phi}$ generated by MRI (see Section 3.5).

3.3. Disk Structure

The time and azimuthally averaged spatial structure inside $30r_{g}$ for the density, flow velocity, radiation energy density and magnetic field lines of the two simulations are shown in Figure 3. Radial profiles of different shell averaged pressure components, the mass accretion rates as well as total optical depth are shown in Figures 4, 5 and 6. The disk is thinner in the run AGN0.07 compared with AGN0.2 due to lower accretion rate and the thickness of the disk is different from standard accretion disk model predictions. In a radiation pressure and electron scattering dominated $\alpha$ disk model in a spherically symmetric gravitational potential, vertical hydrostatic equilibrium implies that the thickness of the disk $H=[\kappa_{s}\dot{M}/(4\pi c)](d\ln\Omega/d\ln r)$ away from the inner boundary. Near the inner boundary $r_{\rm in}$, it is smaller by $1-\left(r_{\rm in}/r\right)^{1/2}$. This is nearly independent of radius (and of $\alpha$), and implies $H/r_{g}\approx 1.1$ (Frank et al., 2002) for $\dot{M}=7\%\dot{M}_{\rm{Edd}}$ and $H/r_{g}\approx 3.0$ for $\dot{M}=20\%\dot{M}_{\rm{Edd}}$. It will be even smaller if we take into account the boundary effect. As shown in Figure 3, the heights of the photosphere clearly increase rapidly with radius in the two runs, and become larger than the values in the standard $\alpha$ disk model beyond $\approx 10r_{g}$. The contrast of radiation energy density inside and above the photosphere is much smaller in AGN0.07 compared with the run AGN0.2 because of the reduced total optical depth as well as different spatial distributions of dissipation. Although the shell averaged radiation pressure is still larger or comparable to the magnetic pressure (Figure 4), the vertical gradient of radiation pressure is actually smaller, particularly in the run AGN0.07. The disk is actually supported by magnetic pressure gradient in this region (section 3.5). There are two magnetic field loops above and below the disk midplane in the run AGN0.07 and they are configured in such a way to have a net radial magnetic field near the midplane. In the run AGN0.2, there are net poloidal magnetic fields threaded through the disk by design. With both magnetic field configurations, we have confirmed that the strong magnetic pressure is dominated by the turbulent component, since the magnetic pressure due to the azimuthally averaged mean magnetic field $\overline{B}^{2}/2$ is smaller than the total magnetic pressure by a factor of $\approx 10$ (Figure 4).

Radial profiles of mass accretion rate averaged over the whole spherical shell reach constant values inside $\approx 15r_{g}$ for AGN0.07 and $\approx 13r_{g}$ for AGN0.2 (Figure 5). We decompose the net mass accretion rate into two parts based on the location of the photosphere shown in Figure 3. At each radius, we integrate $\rho v_{r}r^{2}$ for $\theta$ angles either above or below the photosphere and radial profiles of the two components are shown as red and dashed black lines in Figure 5. For the run AGN0.2, the majority of accretion happens in the optically thick part of the disk until $r<7r_{g}$. While for the run AGN0.07, more than half of the mass is accreted in the optically thin surface of the disk over the entire radial range where the disk has reached steady state. The spatial distribution of mass inflow is consistent with the location of stress for angular momentum transport (section 3.4). Significant surface accretion has also been found in non-radiative ideal MHD simulations of accretion disks with a net poloidal magnetic field (Zhu & Stone, 2018). In fact, the resulting magnetic field configuration from the run AGN0.2 shown in Figure 3 is very similar to what Zhu & Stone (2018) found (for example, Figure 25 of that paper). The poloidal magnetic fields are advected inwards above the disk midplane. Since we can determine the thermal properties of the disk self-consistently, for the accretion rate of $26\%\dot{M}_{\rm{Edd}}$, accretion still happens in the optically thick part of the disk with little accretion in the corona region. When the accretion rate drops to $7\%$ with a reduced optical depth, mass is primarily accreted in the optically thin region. However, for these AGN disks, the mechanism for angular momentum transport responsible for the corona accretion is due to radiation instead of magnetic field as found in non-radiative ideal MHD simulations (section 3.4).

Radial profiles of time averaged Rosseland mean optical depth $\tau_{R}=\kappa_{R}\Sigma/2$ for the two runs are shown in Figure 6. Here the Rosseland mean opacity $\kappa_{R}\equiv\kappa_{a}+\kappa_{s}$ is the sum of absorption $\kappa_{a}$ and scattering $\kappa_{s}$ opacity. Optical depths for effective absorption $\tau_{\rm eff}=\sqrt{\kappa_{s}\kappa_{a}}\Sigma/2$ are also shown in the same Figure. The Rosseland mean optical depth is only 2 at $6r_{g}$ for AGN0.07 and it increases to $100$ at $14r_{g}$, while the effective absorption optical depth drops below $1$ near the ISCO. When the accretion rate is increased to $26\%\dot{M}_{\rm{Edd}}$, the optical depth is increased by a factor of $2-3$ and the entire disk is optically thick even for effective absorption. The optical depth in the disk is smaller than what the standard $\alpha$ disk model predicts by a factor of $\approx 10$ for the same mass accretion rate due to increased disk scale height and inflow velocity because of strong magnetic pressure support (Section 3.5).

3.4. Stresses for Angular Momentum Transport

In accretion disks, transport of angular momentum can be provided by the sum of Maxwell, Reynolds and radiation (the off-diagonal components of the radiation pressure) stresses. Since radiation pressure is so significant in AGN accretion disks, radiation stress can potentially play an important role as an effective viscosity (Loeb & Laor 1992). In the optically thick regime, radiation viscosity is proportional to $E_{r}/(c\rho\kappa_{R})$ multiplied by the shear rate of the velocity field (Weinberg 1971; Mihalas & Mihalas 1984; Kaufman & Blaes 2016). If we only consider the shear rate for a Keplerian disk near the disk midplane, the ratio between the radiation stress and radiation pressure at radius $r$ is $\sim\mathcal{O}\left(\sqrt{r_{g}/r}/(r\rho\kappa_{R})\right)$. For a standard thin accretion disk model where $r\rho\kappa_{R}\gtrapprox 100$, the radiation stress will be much smaller than the typical Maxwell and Reynolds stresses produced by MRI turbulence. In other words, radiation viscosity is small because the mean free path of the photons is too short. This is also true for super-Eddington accretion disks around AGNs due to large optical depth (Jiang et al., 2017a), even though radiation pressure is significantly larger than the gas and magnetic pressure there. However, this is not the case for these sub-Eddington simulations.

Radial profiles of the shell averaged Maxwell ($\alpha_{m}$), Reynolds ($\alpha_{h}$) and radiation stresses (only the $r-\phi$ component remains, $\alpha_{\rm rad}$) scaled with the shell averaged total pressure for the two runs are shown in Figure 7. The stresses are calculated in the same way as in Jiang et al. (2017a). For both runs, we have also calculated the Maxwell stress due to the azimuthally averaged mean $\langle B_{r}\rangle$ and $\langle B_{\phi}\rangle$, which are shown as the dashed red lines in Figure 7. It is smaller than the total Maxwell stress by a factor of 6 in AGN0.2 and a factor of 2 in AGN0.07, despite the fact that the latter has shearing net radial magnetic field in the midplane. The total stress in the run AGN0.2 is dominated by the Maxwell stress with the Reynolds and radiation stresses smaller by a factor of $\approx 5$, which is consistent with the results reported by previous both local and global simulations of MRI turbulence with net vertical magnetic fields (Bai & Stone, 2013; Zhu & Stone, 2018). However, for the run AGN0.07, the radiation stress is larger than the others by more than a factor of $2$ for all the radial range $r>5r_{g}$, while the Maxwell and Reynolds stresses vary from $2\%$ to $10\%$ of the total pressure, which are typical values found by MRI turbulence without net vertical magnetic fields (Hawley et al., 1995; Miller & Stone, 2000; Hirose et al., 2006; Davis et al., 2010; Jiang et al., 2013a, 2014b).

Radiation stress and Maxwell stress dominate different locations of the disk in the run AGN0.07 as they have different spatial distributions. The time and azimuthally averaged spatial distributions of the $r-\phi$ and $\theta-\phi$ components of radiation stress $P_{r}^{r\phi}$ and $P_{r}^{\theta\phi}$ for the run AGN0.07 are shown in the first two panels of Figure 8. For comparison, the Maxwell stress $S_{m}$ and Reynolds stress $S_{h}$ are shown in the third and fourth panels. Inside the photosphere of the disk, the radiation stress is actually much smaller than the Maxwell and Reynolds stresses. This is consistent with our estimate that radiation stress is proportional to $1/(r\rho\kappa_{R})$, while the Maxwell and Reynolds stresses roughly follow the distribution of mass. However, once near the photosphere, because of the increased photon mean free path as well the large velocity gradient, significant radiation stress shows up while Maxwell and Reynolds stresses drop significantly. This suggests that in the run AGN0.07, angular momentum is transported by the Maxwell and Reynolds stresses in the optically thick part of the disk. But once near the photosphere, angular momentum is transported by the radiation stress. Since total radiation stress dominates, more accretion happens in the optically thin region, which is consistent with the mass accretion rates shown in Figure 5. This is not the case for the run AGN0.2, where most accretion still happens in the optically thick region until $r\lessapprox 8r_{g}$. This also causes a stronger corona component at large radii for the run AGN0.07 (Section 3.6), since dissipation associated with the accretion in the optically thin region can heat up the gas easily without cooling efficiently.

The amount of radiation stress we observe in the simulations can also be compared with analytic expectations for radiation viscosity. When the optical depth per cell is large, the $r-\phi$ component of the co-moving frame radiation viscosity can be calculated analytically as (Masaki, 1971; Kaufman & Blaes, 2016)

[TABLE]

where the $r-\phi$ shear rate in spherical polar coordinates is $D_{r\phi}=r\partial\left(v_{\phi}/r\right)/\partial r$+$1/\left(r\sin\theta\right)\partial v_{r}/\partial\phi$. However, when cells become optically thin, this formula needs to be modified to account for the nonzero mean free path of photons. Under optically thin conditions, the viscous stress tensor associated with a shear velocity difference across some given length scale should be proportional to the optical depth across that scale times the velocity difference (Socrates et al., 2004; Kaufman & Blaes, 2016). In order to ensure continuous behavior between the optically thin and thick regimes, we therefore write the stress tensor as

[TABLE]

Here $D^{\prime}_{r\phi}\equiv r\Delta(v_{\phi}/r)+\Delta v_{r}$ is just the difference of shear rate across distances with optical depth $\tau$ along radial (for $v_{\phi}$) and azimuthal (for $v_{r}$) directions. and $k(\tau)=\tau^{-1}$ for $\tau>1$ and $\tau$ for $\tau<1$. In practice, we first determine the optical depth per cell at each location. If that value is larger than one, we calculate the difference of shear rate between neighboring cells. If that value is smaller than 1, we extend the distance until $\tau=1$ is reached and calculate the difference of shear rate there. The radiation stress tensor calculated based on the above formula is transformed back to the lab frame and then averaged azimuthally to compare with the $r-\phi$ radiation stress returned by the simulations directly. Comparisons for two different radii are shown in Figure 9. The numerically calculated radiation stress nicely follow the prediction for radiation viscosity in both optically thick and optically thin regimes.

3.5. Vertical Structure of the Disk

Time and azimuthally averaged poloidal profiles of various quantities at $10r_{g}$ for the two runs are shown in Figure 11. Density drops faster with height in the run AGN0.07 due to a smaller pressure scale height. Shapes of density profiles are also more centrally peaked compared with gas or radiation pressure dominated disks found by both local shearing box and global simulations (Turner, 2004; Hirose et al., 2006, 2009; Jiang et al., 2016, 2017a). Although the shell averaged radiation pressure is comparable to the shell averaged magnetic pressure in the run AGN0.07 (Figure 4), the radiation pressure is relatively flat with height and the whole disk is supported by the magnetic pressure gradient. In fact there is a small enhancement of $P_{r}$ near the photosphere (see also Figure 3) due to the significant dissipation there caused by radiation viscosity. In the run AGN0.2, radiation pressure partially supports the disk near the midplane and the disk becomes completely magnetic pressure supported around $10^{\circ}$ away from the midplane. This also causes the density to drop more slowly with height around the same location as shown in the right panel of Figure 11. Gas pressure is always smaller than the other pressure components by more than a factor of 1000.

For the two magnetic field configurations used in the simulations, besides the Maxwell stress from the turbulence, there are always significant azimuthally averaged mean $\langle B_{r}\rangle$ and $\langle B_{\phi}\rangle$, although the product of these components never becomes the dominant stress (Figure 7). This is different from the magnetic pressure supported disk as found by Gaburov et al. (2012), where the Maxwell stress is primarily due to $-\langle B_{r}\rangle\langle B_{\phi}\rangle$. They also have different vertical distributions compared with the turbulent stress as shown in the third panels of Figure 11. The Maxwell stress generated by the turbulence (the dashed red lines) shows double peaks away from the midplane, which is consistent with previous simulations (Blaes et al., 2011; Jiang et al., 2016). The stress due to the azimuthally averaged mean magnetic field (the difference between the solid black lines and the dashed red lines) is peaked at the midplane and drops quickly with height.

Inside the photosphere near the midplane, the rotation speed of the disk is pretty close to the Keplerian value with negligible radial inflow speed (the fourth and bottom panels of Figure 11). Inside the corona region, the rotation speed drops to only $40\%$ of the Keplerian value and significant radial velocity is present.

Despite the fact that the disk is supported by magnetic pressure, the energy dissipated in the disk is transported away vertically by both the radiation flux and Poynting flux near the midplane. Vertical profiles of the poloidal components of the energy fluxes at $10r_{g}$ for the run AGN0.07 are shown in Figure 12. They are scaled with the critical energy flux $F_{c}\equiv cGM_{\rm{BH}}|\cos(\theta)|/\left(\kappa_{\rm es}(r-2r_{g})^{2}\right)$, which is the radiation flux if the vertical component of gravity is completely balanced by the radiation force with the electron scattering opacity $\kappa_{\rm es}$. The radiation flux is less than $\approx 5\%$ of the critical value inside the disk, which confirms that radiation pressure plays a negligible role to support the disk. But the radiation flux is comparable to the Poynting flux within $\approx 3^{\circ}$ away from the midplane and it completely dominates the energy transport beyond that. The sign of Poynting flux also suggests that magnetic energy amplified in the disk is transported away from the midplane likely due to magnetic buoyancy (Blaes et al., 2011; Begelman et al., 2015), although this simulation does not show any butterfly diagram (Figure 2). The existence of non-zero azimuthally averaged radial magnetic field $\langle B_{r}\rangle$ in localized region of the disk will also keep generating $B_{\phi}$ due to the differential rotation of the disk, which is balanced by the buoyant escape of the field. Since the shell averaged $B_{r}$ is zero, $B_{\phi}$ generated at different locations in the disk will reconnect and cause dissipation.

3.6. Properties of the Coronae

Local shearing box radiation MHD simulations (Jiang et al., 2014b) have found that coronae with high gas temperature are only formed above the disk when the surface density is low enough so that a significant fraction of the total dissipation happens in optically thin regions. For the surface density at $30r_{g}$ adopted by Jiang et al. (2014b) based on the standard thin disk model, only $3.4\%$ of the dissipation happens in the corona region. This fraction is increased to more than $50\%$ inside $15r_{g}$ for the run AGN0.07 with significantly reduced surface density in the magnetic pressure supported disk. This results in a high gas temperature of $\approx 10^{8}$ K covering all the radial range $\lesssim 15r_{g}$ with height $20r_{g}$ away from the disk midplane as shown in the left panel of Figure 10. When the accretion rate is increased to $26\%\dot{M}_{\rm{Edd}}$ in the run AGN0.2, the disk becomes thicker and the fraction of accretion as well as the associated dissipation in the optically thin region are also reduced. Coronae with high gas temperature ($\gtrapprox 10^{8}$ K) only show up inside $\lesssim 10r_{g}$ when the surface density drops to a value such that $\kappa_{R}\Sigma/2\lesssim 100$. Inside this region, the corona is also more vertically extended covering more than $40r_{g}$ away from the disk midplane.

Coronae produced by these simulations are pretty consistent with many observational inferred properties of AGN corona. The simulations only produce high temperature corona (gas temperature $>10^{8}$ K) inside $\approx 10r_{g}$ and the temperature decreases with increasing distance from the black hole. This is because surface density of the disk increases with radius and the fraction of dissipation in the optically thin region drops, which is necessary to maintain the high temperature corona. This agrees with the constrains from both reflection and microlensing modeling of AGN coronae (Chartas et al., 2009; Dai et al., 2010; Reis & Miller, 2013; Jiménez-Vicente et al., 2015), which suggest that they are compact with roughly the same size as we find.

Surface density increases for the same region of the disk when the accretion rate increases from AGN0.07 to AGN0.2. This is in contrast to the standard thin disk model, for which the surface density decreases linearly with increasing mass accretion rate in the radiation pressure dominated regime. When the surface density increases, the fraction of dissipation in the optically thin regime is reduced and the radial regime where the corona is formed is also smaller. This will result in a smaller ratio between X-ray luminosity from the corona and luminosity of the thermal emission from the disk, which is in excellent agreement with the observational fact that the spectral index of AGNs gets harder as luminosity decreases (Steffen et al., 2006; Just et al., 2007; Ruan et al., 2019).

4. Discussion

It is interesting to compare the structures of the disk produced by the simulations with proposed models of magnetic pressure dominated disks in the literature. Pariev et al. (2003) modified the standard $\alpha$ disk model by including a magnetic pressure component in the disk, which was assumed to be larger than the thermal pressure at all radii. Magnetic pressure is assumed to be a constant at each radius and decreases with increasing radius following a power law, the slope of which is a free parameter. The strength of the magnetic field is also a free parameter in this model. The disk structures in this model are very similar to the original $\alpha$ disk model, which predicts that the surface density is proportional to $1/\dot{M}$. This is clearly in contrast to what we find from the simulations. Begelman & Pringle (2007) proposed that the toroidal magnetic field should saturate to the level so that the associated Alfvén velocity is the geometric mean of gas isothermal sound speed $c_{s}$ and Keplerian speed $V_{k}$ . Under this assumption, the disk surface density is found to be proportional to $\dot{M}^{7/9}$, which has the same trend as what we find, although their dependence on $\dot{M}$ is weaker. The Alfvén velocity in our simulations is found to be larger than $\sqrt{c_{s}V_{k}}$ by a factor of $2-3$ in the midplane and much smaller than that near and above the photosphere. Our simulation results can be used to constrain the assumptions and improve these analytical models.

The relatively uniform spatial distribution of $E_{r}$ in the run AGN0.07 is caused by the significant dissipation near the photosphere. In fact, $E_{r}$ near the photosphere is slightly larger than $E_{r}$ at the midplane as shown in the left panel of Figure 11. This is also consistent with the rapid increase of radiation flux at these locations as shown in Figure 12. For comparison, when more dissipation happens in the optically thick part of the disk as in the run AGN0.2, radiation energy density decreases with increasing height monotonically. This is also the reason why radiation viscosity can be significant in the run AGN0.07. Since radiation viscosity is proportional to $E_{r}/\tau$ in the optically thick regime, if $E_{r}$ at the photosphere is smaller than $E_{r}$ in the midplane by a factor of $1/\tau$ as in the standard thin disk model, radiation viscosity will not be significant near the photosphere even mean free path is larger there.

These simulations adopt the gray approximation for radiation transport and Compton scattering is included based on the effective radiation temperature to roughly estimate the energy exchange between photons and gas. Therefore, we cannot produce the expected spectra from these simulations directly. This approximation will also cause underestimate of gas temperature in the corona region. Monte Carlo technique has been developed to post process the simulation data and generate spectra from the disk. Preliminary results show that the coronae in these simulations are indeed able to produce significant X-ray emission with X-ray luminosity consistent with the total dissipation we have in these region. Detailed spectrum properties will be presented in a separate publication. Since corona only exists inside $\approx 10r_{g}$, for which general relativity may play an important role, simulations with fully self-consistent radiation transport in general relativity will be carried out in the near future.

The total luminosity from the inner $10r_{g}$ of the disk is only $1\%L_{\rm{Edd}}$ with $\dot{M}_{\rm{Edd}}=7\%\dot{M}_{\rm{Edd}}$ for the run AGN0.07, which means the radiative efficiency is only $1.4\%$. The radiative efficiency is slightly increased to $\approx 3.5\%$ for the run AGN0.07. The efficiency from these disks is significantly smaller than the prediction of the standard thin disks for the same radial range. One reason for the lower efficiency is the significant dissipation in the optically thin region, which cannot convert to thermal radiation efficiently. Another reason for the lower radiative efficiency in the run AGN0.07 is the stress $-\,\langle\langle B_{r}\,\rangle\rangle\,\langle\langle B_{\phi}\,\rangle\rangle$ near the disk midplane (Figure 11). This component of stress can cause mass accretion without dissipation.

Both of the simulations in this paper established long-lived thermal equilibria between turbulent dissipation of gravitational binding energy in the flow and radiative cooling for more than $10$ thermal time scales. This is despite the fact that radiation dominates the thermal pressure, a situation which generally leads to thermal instability in $\alpha$-disk models (Shakura & Sunyaev, 1976), at least when electron scattering dominates the Rosseland mean opacity. Iron opacity can be significant in disks around supermassive black holes, and has been shown to stabilize the thermal instability in shearing box simulations of MRI turbulence (Jiang et al., 2016). However, while we have included this source of opacity in the simulations here, it actually does not play a significant role in the inner regions near the black hole. Instead, it is likely that the fact that these simulations are supported by magnetic, not thermal pressure, leads to their stability (Begelman & Pringle, 2007). Indeed, Sa̧dowski (2016) found that global simulations of geometrically thin disks are stabilized if the magnetic pressure exceeds half the total pressure, a criterion which we more than satisfy in our simulations here.

Light curves of a variety of accretion-powered sources in astrophysics, including active galactic nuclei, show a linear relationship between rms variability amplitude and flux (see e.g. Uttley et al. 2005 and references therein). Thereis a hint of this in the AGN0.2 light curve shown in Figure 1, where the variability amplitude is clearly larger at higher luminosities than at lower luminosities. Unfortunately, the simulation was not run long enough and the light curve is too short to produce a clean rms-flux relation. On the other hand, there is no indication of any relationship between variability amplitude and flux in AGN0.07. Hogg & Reynolds (2016) have suggested that the rms-flux relation is related to the presence of the characteristic butterfly diagram of quasi-periodic azimuthal field reversals. It is therefore noteworthy that these are present only in AGN0.2, and not in AGN0.07 because of its much stonger magnetic pressure support. The presence of an rms-flux relation may be a way of observationally ruling out magnetic fields which are so strong that the butterfly dynamo is suppressed, but longer simulations will need to be run to confirm this.

Note that the shape and magnetic field in the initial tori we use for the two simulations AGN0.07 and AGN0.2 are very similar to the ones adopted by Jiang et al. (2017b). We only change the density of the torus to achieve different mass accretion rates, while the ratio between the initial magnetic pressure and radiation pressure is kept fixed. However, the accretion disks with super-Eddington accretion rates shown in Jiang et al. (2017b) do not evolve to a magnetic pressure dominated regime. When the surface density as well as the optical depth are increased, radiation pressure becomes more important in supporting the disk. Since the amount of magnetic field that in the disk is determined by a balance between buoyant escape and amplification by the MRI, this suggests that stronger radiation support in the disk may either increase the magnetic buoyancy or reduce the saturation amplitude of MRI (Jiang et al., 2013b). We intend to investigate this further in future work.

5. Summary

We have successfully simulated two accretion disks around a $5\times 10^{8}M_{\rm{BH}}$ black hole with mass accretion rates reaching $0.07\dot{M}_{\rm{Edd}}$ and $0.2\dot{M}_{\rm{Edd}}$ up to $\sim 15r_{g}$. The disks do not show any sign of thermal instability over many thermal timescales. The structure of the disk differs markedly from the standard thin disk model (Shakura & Sunyaev, 1973) for the same mass accretion rate in the following ways.

• •

The disk is supported vertically by magnetic pressure rather than thermal pressure for these accretion rates.

• •

The surface density and total optical depth are reduced by a factor of $\approx 10$. The disk scale height increases with radius significantly instead of being constant.

• •

A significant fraction of dissipation as well as the associated mass accretion (more than $50\%$ for $\dot{M}=0.07\dot{M}_{\rm{Edd}}$) happen away from the disk midplane in the optically thin region, resulting in a high temperature corona inside $10r_{g}$. The fraction of dissipation in the corona region decreases with increasing mass accretion rate.

• •

The transport of angular momentum is due to Maxwell and Reynolds stress inside the photosphere. But near the photosphere, radiation viscosity is the dominant mechanism for angular momentum transport, and this can dominate the vertically averaged stress.

Acknowledgements

We thank Julian Krolik for helpful comments on an earlier draft. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915, 17-48958 and AST-1715277. S.W.D. is supported by a Sloan Foundation Research Fellowship. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Resources supporting this work were also provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.