Magneto Rotational Instability in Magnetized AGN Tori

TL;DR

This study uses 3D MHD simulations to explore how magnetic fields and MRI influence mass supply and structure in AGN tori, revealing magnetic field reversals and resolution requirements for MRI.

Contribution

It demonstrates the role of MRI in AGN tori dynamics and establishes the necessary numerical resolution to accurately simulate MRI effects.

Findings

Magnetized disks develop cold and warm gas components.

Azimuthal magnetic fields reverse direction quasi-periodically.

MRI wavelength must be resolved with at least 20 grid cells.

Abstract

It is widely believed that in active galactic nuclei (AGNs) a supermassive black hole with an accretion disk is surrounded by an optically and geometrically thick torus at sub-parsec scale. However, it is not clear how is the mass supply toward the central engine caused and how it is related with the internal structures of the tori. The magnetic field in the tori may contribute to the accretion process via the magneto-rotational instability (MRI). Using global three dimensional magnetohydrodynamic (MHD) simulations taking the effects of X-ray heating and radiative cooling into account studied the numerical resolution for azimuthal direction for MRI driving. We found that a strongly magnetized disk consisted of a cold ($< 10^3$ K) and warm ($10^4$ K) gas is developed in about 30 rotational periods. We also found in high resolution model that the mean azimuthal magnetic fields reverse their direction quasi-periodically. We confirmed that the typical wave length of the MRI should be resolved with a least 20 azimuthal grid cells.

Figures

Equations

$$\dot{M}=\frac{L}{\eta c^2}\sim 0.2\left[\frac{M}{10^7{M}_{\odot }}\right]M_{\odot }{\mathrm{yr}}^{-1},$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left[\rho \mathbf{v}\right]=0$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left[\rho \mathbf{vv}+\left(P_{\mathrm{g}}+\frac{B^2}{8\pi }\right)\mathbf{I}-\frac{\mathbf{BB}}{4\pi }\right]=-\rho \nabla \Phi$$

$$\frac{\partial }{\partial t}\left(\frac{P_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}+\frac{1}{2}\rho v^2+\frac{B^2}{8\pi }\right)+\nabla \cdot \left[\left(\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}{P}_{\mathrm{g}}+\frac{1}{2}\rho v^2\right)\mathbf{v}-\mathbf{E}\times \mathbf{B}\right]$$

$$=-\rho \mathbf{v}\cdot \nabla \Phi +\rho L$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \mathbf{E},$$

$$\mathbf{E}=\mathbf{v}\times \mathbf{B}-\eta \nabla \times \mathbf{B},$$

$$\rho L=n\left({\Gamma }_{\mathrm{UV}}+{\Gamma }_{\mathrm{Coulomb}}\right)+n^2\left({\Gamma }_{\mathrm{Compton}}+{\Gamma }_{\mathrm{photoionic}}-\Lambda \right),$$

$${\Gamma }_{\mathrm{Compton}}=8.9\times 10^{-36}\xi \left(T_X-4T\right)[\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^3],$$

$${\Gamma }_{\mathrm{photoionic}}=1.5\times 10^{-21}{\xi }^{1/4}{T}^{1/2}\left(1-T/T_X\right)[\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^3],$$

$${\Gamma }_{\mathrm{Coulomb}}={\eta }_{\mathrm{h}}{H}_{\mathrm{X}}[\mathrm{erg}{\mathrm{s}}^{-1}],$$

$$\xi =\frac{L_{\mathrm{X}}}{n\left(r^2+z^2\right)}\exp (-\tau )[{\mathrm{erg}\mathrm{cm}\mathrm{s}}^{-1}],$$

$$Q_{\phi }=\frac{{\lambda }_{\mathrm{MRI}}}{r\Delta \phi }=\frac{2\pi }{\Delta \phi }\frac{|B_{\phi }|}{v_{\phi }\sqrt{4\pi \rho }},$$

Full text

galaxies: nuclei — methods: numerical — magnetohydrodynamics (MHD)

††journal: Galaxies of the Special Issue The Power of Faraday Tomography

1 Introduction

Accretion processes onto central objects are important for evolution of the astrophysical objects; e.g. the X-ray binary, Gamma ray burst, active galactic nuclei (AGNs). Since the most accreted gases have angular momentum, angular momentum transport is a problem. Turbulent viscosity was suggested by Shukura $\&$ Sunyaev (1973), but driving source of the turbulence was not identified. Balbus $\&$ Hawley (1991) pointed out that the magneto-rotational instability (MRI) can be account for the angular momentum transport in the differentially rotational magnetized disks. Mass accretion and angular momentum transport in magnetized gas disks can be described as follows. Small amplitude perturbation of the magnetic field lines in radial direction are stretched and amplified for the radial direction by MRI growing. Since magnetic field lines are frozen into the ionized gas under the MHD approximation, deforming magnetic fields transport the gas and the angular momentum to central region. A viscosity proportional to the gas pressure was conventionally assumed in the equation of motion, but it is replaced by Maxwell’s stress, $B_{r}B_{\varphi}$.

MRI-driven turbulence amplifies the magnetic field strength in the linear regime, and drives a buoyancy force due to Parker instability (Parker, 1966). Nonlinear evolution of MRI with the Parker instability is characterized by a quasi-periodic reversal of the direction of azimuthal field in spacetime diagrams (e.g., Beckwith et al. (2011); Machida et al. (2013); Parkin et al. (2013); Hogg $\&$ Reynolds (2016)). The typical growing timescale of MRI is the rotational timescale determined by a balance of radial direction between the gravity, centrifugal force and magnetic tension, and the period of the quasi-periodic reversal is about 10-20 rotational periods. Therefore, a long-term calculation beyond 10 rotational periods are required in order to study the time evolution of nonlinear evolution of MRI and the mass transfer due to the MHD turbulence.

Here we focus on sub-pc structures of the magnetized gas around a supermassive black hole (SMBH). The accretion rate in the AGNs is related to their luminosities as

[TABLE]

where $L$ is the bolometric luminosity, $c$ is the light speed, and $\eta$ is efficiency of the energy conversion, respectively. The SMBHs are surrounded by optically and geometrically thick tori (see e.g., Antonucci (1993); Urry $\&$ Padovani (1995)), and the accreted material toward SMBH could be originated in the tori. Magnetic fields penetrating the torus were reported by the mid-infrared spectro-polarimetry of NGC1068 and NGC4151 (Lopez-Rodriguez et al., 2016, 2018). Therefore, for sake of understanding the torus accretion, MRI is one of the mechanisms of the angular momentum transport. Dorodnitsyn $\&$ Kallman (2017) and Chan $\&$ Krolik (2017) carried out the simulations of magnetized AGN tori. However the nonlinear MRI insufficiently have understood. It is not clear whether MRI drives with cold gas in AGN tori.

We have studied the MRI in the cold gas of AGN torus using global three-dimensional MHD simulations including the X-ray heating and the radiative cooling. We investigate how the azimuthal resolution depends on driving MRI in the cold dense gas using the high-order numerical scheme.

2 Methods

We used CANS+ code (Matsumoto et al., 2016), which is implemented with the HLLD method (Miyoshi $\&$ Kusano, 2005) and the hyperbolic divergence cleaning method (Dedner et al., 2002). and the 5-th order spatial accuracy is achieved by the monotonicity preserving method (MP5 of Suresh $\&$ Huynh (1997)). The high order interpolation requires to reduce the numerical diffusion of magnetic fields for a long-term calculation. A cylindrical coordinate $(r,\varphi,z)$ is used in the computational domain, $0<r~{}[{\rm pc}]~{}<11$, $0<\varphi<2\pi$, and $|z|~{}[{\rm pc}]>3$. The numbers of grid points are $(N_{r},N_{z})=(256,512)$, and we chose $N_{\varphi}$ = 128 or 512. For the outer region ($|z|>2.8$[pc], $r>10$ [pc]) and the central region ($\sqrt{r^{2}+z^{2}}<0.4$), the outflow boundary conditions and absorbing boundary conditions are used, respectively.

We solve the following MHD equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$, $P_{g}$, $\bm{v}$, $\bm{B}$ is the gas density, pressure, velocity vector, magnetic field, respectively. In order to evaluate the temperature, we assumed the ideal gas with $\gamma_{\rm g}=5/3$. Electric field, $\bm{E}$ is related as the Ohm law, and the magnetic resistivity $\eta$ adopts the the anomalous resistivity model (e.g. Yokoyama $\&$ Shibata (1994); Machida et al. (2013)). This is in effect where the magnetic reconnection occur and regulated to be $\eta<10^{4}(r/1~{}[{\rm pc}])(v/207~{}[{\rm km/s}])$ [cm2/s] The gravitational potential $\Phi$ is $GM/\sqrt{r^{2}+z^{2}}$, where $M$ is the mass of SMBH, $10^{7}M_{\odot}$. We ignore the self-gravity of the gas.

We combined the radiative cooling and heating effects into eq. 4,

[TABLE]

where $n\left(=\rho/m_{\rm H}\right)$ is the number density, and $m_{\rm H}$ is the mass of neutral hydrogen. Cooling function was modeled on Wada et al. (2009) and Wada (2012) with the solar abundances. We take X-ray and UV from the accretion disk into account as heating processes. $\Gamma_{\rm UV}=1.8\times 10^{-25}[{\rm erg~{}}{\rm s}^{-1}]$ is assumed. X-ray heating (Blondin, 1994) due to Compton interaction and photoionization for $10^{4}<T<10^{8}$ are:

[TABLE]

[TABLE]

where $T_{X}=10^{8}$ K is the characteristic temperature of X-ray. The Coulomb heating is,

[TABLE]

$\eta_{\rm h}$ denotes the efficiency (Dalgarno et al., 1999; Meijerink $\&$ Spaans, 2005), which assume to be fixed $0.2$, and $H_{\rm X}$ is X-ray energy deposition rate $H_{\rm X}=3.8\times 10^{-25}\xi{\rm~{}erg~{}s}^{-1}$. The X-ray luminosity in the nucleus is parameterized by ionization parameter,

[TABLE]

where $L_{\rm X}$ is X-ray luminosity and $\tau$ is the optical depth, respectively. We assumed $\xi=1.0$ for simplification.

Initial condition is expressed as the superposition of dynamical equilibrium solution with the isothermal spherical symmetry and the weakly magnetized hot torus (see, Okada et al. (1989); Machida et al. (2013)). We adopted that the torus center is at 1 [pc] and rotation velocity is $v_{\varphi}=207(r/1~{}[{\rm pc}])^{-0.65}$ [km/s] for the torus and $v_{\varphi}=0$ for otherwise. The initial magnetic field in the torus is $(B_{r},~{}B_{\varphi},~{}B_{z})=(0,~{}P_{g}/\beta,0)$, and plasma beta is 100 defined as $\beta=2P_{\rm g}/\left(B_{r}^{2}+B_{\varphi}^{2}+B_{z}^{2}\right)$. The simulations used the normalized unit, i.e., $l_{0}=1$ [pc] for the length, $v_{0}=\sqrt{GM/l_{0}}=207$ [km/s] for the velocity, $n_{0}=10^{2}$ [cm*-3*] for the number density, and $B_{0}=\sqrt{4\pi m_{\rm H}n_{0}v_{0}^{2}}=3.3$ [mG] for the magnetic field, respectively.

We initially evolve the model adiabatically, and once the MRI-driven turbulence is fully developed and becomes quasi-stable at about $20$ rotational periods at $r=1$ [pc], the cooling and heating terms are taken into account.

3 Results

Fig.1 from left to right shows time evolution in the high resolution model ($N_{\varphi}=512$): (1) initial state, (2) the adiabatic phase where the MHD turbulence is developed, and (3) the cooling/heating phase, respectively. Initial torus of (1) constitutes the high temperature ($10^{4}<T<10^{5.3}$) gases. (2) shows that the gases are heated by Joule heating and spread out the radial and vertical direction from initial torus. The radial spread of gas is a result of the angular momentum transport with the deformation of magnetic field lines attracted by MRI. The vertical spread is caused by the magnetic field to the vertical direction. After the heating and cooling tern on (Fig.1(3)), a geometrically thin disk consisted of cold gas ($<10^{3}$ K), which is surrounded by a warm gas halo ($\sim 10^{4-5}$ K) is formed.

The magnetic field of the $rz$ plane slice is shown in Fig.2: (a) $B_{\varphi}$ and (b) plasma beta. As shown in Fig.2(a2), a turbulent structure is developed due to MRI, It is also notable that the turbulence is consisted of opposite directions of $B_{\varphi}$ represented by blue and red regions. The plasma beta decreases to order unity from the initial value ($\beta$ =100). After cooling and heating are taken account, the spatial sign reversal of $B_{\varphi}$ is dissipated and formed a stripe-like structure on the vertical direction (Fig.2(a3)). As shown by the plasma beta (Fig.2(b3)), the magnetic field is enhanced to $\beta\sim 0.1$ by compression of the gas in a few rotation periods. However, re-amplification to $\beta\sim 0.01$ occur by no compression effects. We will recall this to show in Fig.5.

In the cold, thin disk seen in Fig.1(3) and Fig. 2(a3), the turbulent structures are not apparent. We here investigate whether this is due to lack of the numerical resolution for the azimuthal direction to resolve MRI. Hawley et al. (2013) suggested that Q-value, which is the ratio of grid size to the characteristics MRI wavelength,

[TABLE]

should be large enough to resolve the MRI, e.g. $Q_{\varphi}\geq 20$. Fig.3 shows the $Q_{\varphi}$ distribution in the two models with different spatial resolutions. Low resolution model ($N_{\varphi}=128$) in the top of Fig.3 holds the region of $Q_{\varphi}<20$ regardless of whether the radiative cooling and heating are effective or not. On the other hands, $Q_{\varphi}>20$ in most regions in the high resolution model ($N_{\varphi}=512$). However, one should note that the MRI may be still not well resolved in the cold disk ($1<r~{}[{\rm pc}]<4$), where $Q_{\varphi}\sim 4-15$.

Time evolution of azimuthally averaged $B_{\varphi}$ at $r=1$ [pc] is shown in Fig.4. Red and blue colors represent that the direction of $B_{\varphi}$ is opposite. Before the heating and cooling tern on ($t<0.60$ [Myr]), the direction of $B_{\varphi}$ is varied from positive to negative inside the torus. This implies that the inner magnetic field escapes buoyantly from the torus in the vertical direction due to Parker instability. Although the characteristic timescale growing linear MRI is rotational period, the timescale of the direction reversal appears about 10 rotational periods, $T_{\rm cycle}\sim 2.98\times 10^{-2}$ [Myr] at $r=1$ [pc]. This implies that the quasi-periodic reversal is dominated by the non-linear growth of MRI, as seen in the simulations of the galactic or accretion disks (e.g., Machida et al. (2013); O’Nell et al. (2011); Flock et al. (2012); Parkin et al. (2013); Hogg $\&$ Reynolds (2016)). When the heating and cooling are taken into account ($t>0.60$ [Myr]), the difference between the two models is evident. In the high resolution model, the direction of $B_{\varphi}$ at the mid-plane is reversed, as shown by the change of color from blue to red, at $t\sim 2.0$ [Myr], This shows that the magnetic field escapes from the disk plane vertically. However, the low resolution model does not begin to occur the reversal. The difference of the magnetic field structure is responsible for the azimuthal resolution.

We reveal the relation between the magnetic field structure and the spatially spread due to the gas motion. We show the 2D histograms of the total magnetic field strength and number density in Fig.5. The color contour denotes the mass occupied in the cell, $\int M(\rho,B_{\rm total})~{}d\rho dB_{\rm total}$. Left panels are the snapshots of $t=0.60$ [Myr] before heating and cooling are included. Most of the gas is in the regime where the total magnetic field strength is proportional to the number density on the slope, $B_{\rm total}\propto n$. Assuming the conservation of mass and magnetic flux in the torus, this relation means that magnetized torus is spread-out in the $rz$-plane, as shown in Fig.1(1) and Fig.2(1b). Compression of the cold gases under the cooling effects is shown in middle panels ($t=1.10$ [Myr]). The maximum number density increases up to $>10^{4}$ [cm*-3*], and the field strength is amplified. In the right panels ($t=1.53$ [Myr]), the low resolution model ($N_{\varphi}=128$) indicates that the field strength dose not depend on the number density, $B_{\rm total}\propto n^{0}$. This means that the magnetic field structure leaves unchanged. This is agreement with the left panel of Fig.4. On the other hand, the high resolution model ($N_{\varphi}=512$) are formed the field amplification with the relation, $B_{\rm total}\propto n^{1/2}$. The right of maximum field strength is stronger than middle by a factor of $\sim 2$ The higher the resolution, the stronger the magnetic field is generated in the dense region ($1<n<4$).

We measure the inflow rate evaluated with the mass flux passing through the cylindrical surface at each radius (Fig.6). Both model show that the inflow rate at $r=1$ [pc] increases, starting from $1.00$ [Myr] ($N_{\varphi}$=128) or $1.20$ [Myr] ($N_{\varphi}=512$). The inflow rate is larger with the higher spatial resolution. After $t=1.50$ [Myr], the inflow rate at each radius become decreasing in the low resolution model. On the contrary, the high resolution model only shows that the larger the radius, the later inflow rate changes. These results suggest that the angular momentum transport due to the MRI-driven turbulence, and the resultant mass accretion toward the center, are not well resolved for the model with $N_{\varphi}=128$.

4 Conclusions

We performed three-dimensional MHD simulations including heating and cooling effects in the gas around an AGN. We especially focused on the development of the MRI in a cold gas disk. The MRI-driven turbulence is developed and it makes the disk geometrically thick. We found that the azimuthal numerical resolution affects significantly driving MRI. For example, quasi-periodic reversal of the mean azimuthal magnetic field does not occur in the low-resolution model, where the Q-value (Eq. 12) is less than about 20 in the whole region. We confirmed that $Q_{\varphi}$ should be larger than 20, in order to resolve the MRI and buoyantly escape vertically due to Parker instability.

We thanks the anonymous reviewers, Ryoji Matsumoto, and Mami Machida for helpful comment. This work was supported by JSPS KAKENHI Grant Number 16H03959. Numerical computations were carried out on Cray XC30 and XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. Visualization was performed using Python at https://www.python.org/