Multiphase Circumnuclear Gas in Low $\beta$ Disk: Turbulence and Magnetic Field Reversals

We studied the magnetic field structures and dynamics of magnetized multi-phase gas on pc scales around supermassive black holes by using global 3D magnetohydrodynamics (MHD) simulations. We considered the effect of radiative cooling and X-ray heating due to active galactic nuclei (AGNs). The gas disk consists of a multi-phase gas with: (1) cold ($\leq 10^3$ K) and thin, (2) warm ($\sim 10^4$ K) and thick components with a wide range of number densities. The turbulent magnetic energy at maximum is comparable to the thermal and turbulent kinetic energies in the turbulent motion. We confirmed that the turbulent velocity of the warm gas in the ambient cold gas is caused by magnetoconvective instability. The turbulent magnetic field due to magnetorotational instability (MRI) is developed in the disk, but the mean toroidal magnetic field dominates and supports in a quasi-steady state, where the plasma-$\beta$, the ratio between gas pressure and magnetic pressure, is low ($\beta<1$). As often seen in adiabatic MHD simulations of rotating disks, the direction of the mean toroidal field periodically reverses with time even in multi-phase gas structures. The direction reversal is caused by vertically escaping the magnetic flux from the disk and by the combination of the MRI and the Parker instability.


INTRODUCTION
Magnetic field plays a crucial role in gas dynamics and multi-phase gas structures in the central regions of galaxies where various dynamical structures, such as outflows, jets, and turbulent motions of interstellar medium are present on a wide dynamic range from the accretion disk to kpc scales. It is suggested that the magnetic field in the central region of our Galaxy is a few tens µG to mG (see, e.g., Ferrière 2009;Crocker et al. 2010;Han 2017). In particular, Hsieh et al. (2018) estimated the milligauss toroidal field and plasma-β, β = 0.01 − 1, where β is the ratio between thermal pressure and magnetic pressure, Nishiyama et al. (2010Nishiyama et al. ( , 2013 found that the mean toroidal field extends to scale-heights with the galactic latitude |b| < 0.4 • deduced from the Fe 6.7 keV line emission. The circumnuclear magnetic field is also observed in some galaxies. In the proto-typical type-2 Seyfert, NGC 1068, the polarimetry in the infrared suggested that the magnetic field is dominated by a toroidal field with a few ten mG and β ∼ 0.15 (Lopez-Rodriguez et al. 2015). Using the water vapor masers in NGC 4258 (Modjaz et al. 2005), the upper limit of the toroidal magnetic field ∼ 100 mG is inferred.
Global magnetohydrodynamic (MHD) simulations of a rotating gas disk suggested that a turbulent field is developed due to magnetorotational instability (MRI; Balbus & Hawley 1991). Using adiabatic MHD simulations for the Galactic center, Machida et al. (2009) and Suzuki et al. (2015) reported that magnetic turbulence is driven by MRI. Turbulence contributes to the angular momentum and mass transport in the disk. In the nonlinear phase of MRI, it is characterized by a quasi-periodic reversal of the direction of the mean toroidal field (e.g., O'Neill et al. 2011;Beckwith et al. 2011;Flock et al. 2012;Parkin & Bicknell 2013;Machida et al. 2013;Hogg & Reynolds 2016). This field reversal is caused by the vertical transport of the mean field buoyantly from the mid-plane. As a result, the amplification of the magnetic turbulence is saturated and, consequently, the field strength and the plasma-β in the mid-plane are limited in the non-linear phase.
When β 1, the time scale of the direction reversal is often observed to be about ten rotational periods. On the other hand, for β < 1, this time scale is longer in the nonlinear regime depending on β. Local 3D shearing box simulations reported that the reversal pattern is irregular for the isothermal stratified gas disk (Bai & Stone 2013;Salvesen et al. 2016b). They also showed that when β 0.4, the field does not show reversal within 150 rotational periods. Using global MHD simulations, Zhu & Stone (2018) and Mishra et al. (2019) studied the MRI in a low β disk. Contrary to local MHD simulations, the plasma-β is larger than unity at the mid-plane, and no field reversal is observed within 50 rotational periods. Fragile & Sadowski (2017) conducted simulations starting from a strong toroidal magnetic field (β = 0.1), and they found that magnetic field dissipates to β ∼ 10 in a steady state after 10 rotational periods. Salvesen et al. (2016a) pointed out that the poloidal field is necessary to form a strongly magnetized disk. Begelman & Pringle (2007) studied the low β disk by compiling the typical unstable condition for MRI, and showed that even if there was no vertical magnetic field, the low-β MRI is limited by the unstable condition β c s /v K , where v K and c s are the Keplerian and sound speeds, respectively. Therefore, for a given v K , the plasma-β driving MRI is determined by the disk temperature, and it is expected that an MRI with a low β can develop in a low temperature gas disk.
In most of the previous global MHD simulations of the magnetic field in a rotating disk, the gas was assumed to be adiabatic; therefore, the multi-phase nature of the magnetized gas was not well studied. However, the circumnuclear gas on pc scales should consist of cold (T 10 3 K), warm (T ∼ 10 4 K), and hot (T 10 5 K) gases (see, e.g., for AGN: Netzer 2015, and for the Galactic center: Liu et al. 2013). Using global 3D simulations and considering radiative cooling and heating effects, Wada et al. (2009) and Wada (2012) showed that the gas inside tens of parsecs from SMBHs becomes multi-phase. They applied their model to the circumnuclear disk in the Circinus galaxy, which is the nearest type-2 Seyfert galaxy, and found that the model is consistent with multi-wavelength observations (Wada et al. 2016), such as CO(3 − 2) and [CI](1 − 0) emission lines Izumi et al. 2018). However, the magnetic field was not taken into account in their models.
There are several numerical studies for the magnetized circumnuclear gas on a pc scale. Chan & Krolik (2017) and Dorodnitsyn & Kallman (2017) studied the evolution of warm gas with β = 1. However, the timescale of their simulations is about 10 rotational periods, which is not long enough for the magnetic field to become a steady-state. The long-term steady-state behavior of the direction reversal of the mean magnetic field was not well studied. Moreover, MRI with low β in a non-linear phase was not studied for the multi-phase gases, especially below 10 4 K. For example, it is not clear whether turbulence is maintained even in the cold gas, and how different are the structures of the magnetic field compared to the adiabatic gas. The structure of the circumnuclear gas with a magnetic field should be important to consider the effect of the radiation feedback from the AGN. In order to clarify these questions, we study the long-term behavior of the strong magnetized gas around an SMBH by taking into account realistic cooling and heating processes. This paper is organized as follows. In §2, we present the basic equation, initial condition, and numerical model. Cooling and heating processes are described in §2.1. Numerical results are shown in §3. Development of the MHD turbulence is show in §3.1. Time evolution of the mean toroidal magnetic field and the physical origin of the turbulence are discussed in § §3.2 and 3.3. The thermal structures of the magnetized multi-phase gas are presented in §3.4. In §4, we discuss the direction reversal with low β ( §4.1) and the radiation pressure ( §4.2). Finally, we summarize the results in §5.

Basic Equations
We study the pc-scale magnetized gas disk using 3D MHD simulations, considering radiative cooling and various heating effects in the cylindrical coordinate (R, ϕ, z). The resistive MHD equations are: where B is the magnetic field, ρ is the gas density, ρ ≡ m H n ∼ m H (n n + n i ) with number densities of neutral and ionized hydrogen. P g is the thermal gas pressure. The gas temperature is adopted as that of an ideal gas with the specific heat ratio γ g = 5/3. The gravitational potential is assumed to be the Newtonian potential, Φ = GM BH /r, where G is the gravitational constant, M BH = 10 7 M ⊙ is mass of the SMBH, and r = √ R 2 + z 2 is the distance from the SMBH. The electric field E obeys the Ohm's law, E = −v × B + η∇ × B. We assume the anomalous resistivity η as modeled by Yokoyama & Shibata (1994), where v d = |∇ × B|/ρ is the electron-ion drift velocity, v c = 10 8 cm s −1 is the critical velocity, and we adopted η 0 = 10 −9 and η max = 10 −6 pc 2 yr −1 , respectively. Note that the results presented below are not sensitive to η 0 and η max (see also, Appendix A). The radiative cooling and heating term ρL in Equation 5 is given as where the cooling function Λ (Figure 1) is taken from Meijerink & Spaans (2005) and Wada et al. (2009). As a major heating source, we consider that X-ray photons come from the accretion disk, and the heating function is Γ X = Γ Coulomb + Γ Compton + Γ photoionic . The Coulomb interaction is given by where η h is the efficiency in fixed 0.2, and H X is the X-ray energy deposition rate, H X = 3.8 × 10 −25 ξ erg s −1 . For the Compton and photoionization interactions, we use a formula given by Blondin (1994), i.e.

Initial Condition and Normalized Unit
We start the simulation from an equilibrium torus. In order to construct the equilibrium torus solutions with a toroidal field, we assume the polytropic relation, where K is the polytropic constant, and the weak toroidal field B ϕ is given by plasma-β, β = 100. To give the density distribution, we assume the rotation velocity to be the radial distribution of angular momentum L(R), where L 0 (≡ R 0 v 0 ) and R 0 adopt the normalized units and a = 0.35 is the power index. The flux surface Ψ(R, z) is then given as where c 2 s = γ g P g /ρ is the square of the sound speed and v 2 A = B 2 ϕ /(4πρ) is the square of the Alfvén speed. We defined the normalized unit as the torus density maximum, Ψ 0 = Ψ(R 0 , 0). For the condition of positive pressure, we can obtain a torus density distribution Ψ(R, z) = Ψ 0 , We assume that the torus is embedded in the isothermal non-rotating and non-magnetized halo, with pressure and density given by The normalized quantities are listed in Table 1. The polytropic indexes K and K h are parameterized by the square of the velocity ratio, ε = Kρ in the normalization, and we set ε = 0.04 in the torus and ε h = 2 in the halo, respectively.

Numerical Methods and Model Setup
We use CANS+ code (Matsumoto et al. 2019), which is implemented using the HLLD solver (Harten-Lax-van Leerdiscontinuitues; Miyoshi & Kusano 2005) with div B cleaning (Dedner et al. 2002) and three-stage total variation diminishing Runge-Kutta time integration (TVDRK). The fifth-order accuracy in space is achieved through the Monotonicity Preserving method (MP5; Suresh & Huynh 1997) to capture small-scale magnetic fluctuations. Basic equations are solved using conservative forms with geometrical and gravitational source terms explicitly included. After TVDRK updating, the cooling and heating source terms are treated in the implicit operator splitting approach.
We assume the outflow boundary condition for the outer boundaries (i.e. R/R 0 = 10 and |z|/R 0 > 2.8). For the azimuthal direction, the periodic boundary condition is assumed. Meshes around the cylindrical axis are sent to the opposite computational domain of the azimuthal direction, which means that fluid can flow across the polar axis (R = 0). In the central region for r/R 0 < 0.4, we impose an absorbing boundary condition, i.e.
where q and q init are the primitive variables of the TVDRK updating state and the initial state, respectivley. Damping function D(r) is modeled as, We start the simulations assuming the adiabatic MHD, then after the magnetic field strength sufficiently developed and the system becomes a quasi-steady state at t ≤ 0.477 Myr (25 rotational periods at R = 1 pc), the cooling and heating terms are considered until t ≤ 3.64 Myr (97 rotational periods). Figure 2 shows three snapshots (t = 0.000, 0.745 and 2.353 Myr) of the gas temperature and the plasma-β on a R-z plane. Figure 2a show the initial conditions. As explained in Section 2.3, the system evolves adiabatically until t = 0.477 Myr. During this period, gas spreads out vertically, and intense magnetic field fluctuations are developed by the MRI and β > 0.6 inside the torus. The magnetic field is stronger near the surface of the torus with β ∼ 0.2. MRI causes not only turbulent motion but also heating due to magnetic reconnection. After cooling and heating are taken account, the structure of the torus changes ( Figure 2c). The plasma-β around the mid-plane becomes smaller with 0.003 < β < 6.0. The torus becomes geometrically thinner, and it consists of two components: cold disk (T 10 3 K) and warm disk (T = 10 3−5 K). The cold, thin disk (R < 2 pc) is supported vertically by the magnetic field as discussed below.  Figure 3 shows the magnetic field structure at t = 0.744 Myr of the adiabatic MHD state. The turbulent field in the torus dominates in the R-z and R-ϕ planes. The toroidal field (B ϕ ) shows the flux bundle of the positive direction (red) around the torus surface at |z| = 1.5 pc. Figure 4 is the same as Figure 3, but cooling and heating are considered (t = 2.353 Myr). The field strength and turbulent fluctuation are markedly different from those in Figure 3. The toroidal field B ϕ dominates the total magnetic field. The two plots of B ϕ show that the mean toroidal field has more coherent structures compared to the that in adiabatic phase ( Figure 3) with opposite directions shown in blue and red in Figure 4. On the other hand, radial and R-z plane vertical magnetic fields are dominated by the turbulent component. We can also see that the patches of B R and B z with opposite directions tend to extend radially and vertically, respectively. For a more quantitative observation, we decomposed the magnetic energy into mean and turbulent components. We measured the mean component as the azimuthal average,

Development of the MHD-turbulence in the Torus
Hence, the turbulent component is derived as, where we notate the mean as¯and the turbulence as δ. Figure 5 shows the evolution of the magnetic energy in each component. It is clear that the structures of the magnetic field change drastically after the cooling and heating. For t < 0.5 Myr, turbulent fields become exponentially stronger than mean fields, while the ratios between them are approximately constant, i.e. δB 2 ϕ /δB 2 R ∼ 10 and δB 2 ϕ /δB 2 z ∼ 22. After cooling/heating at t = 0.744 Myr, the turbulent fields decrease quickly by one order of magnitude, and they survive until t ∼ 3.5 Myr. This turbulent field dissipation is caused by the decrease in the gas temperature (see in Section 4). During this phase, we found that |δB ϕ | > |δB R | > |δB z |, but the ratios between them change.
In the adiabatic MHD state, the amplitude of B ϕ (blue solid line in Figure 5) reaches a quasi-steady state with a periodical cycle of T cycle = 0.163 Myr ∼ 9.1 rotational periods. After cooling/heating effects are included, oscillating amplification of B ϕ continues for 0.744 < t < 1.75 Myr. The quasi-steady state is achieved when B ϕ does not significantly change . However, |B ϕ | shows quasi-periodic oscillations on a long time scale beyond reaching steady state. This amplitude and period become larger and longer than those of the adiabatic MHD state. Around the maximum of |B ϕ | (i.e. t = 1.75, 3.08 Myr), |δB| and |B R | are maximized, and around the minimum of that (i.e. t = 2.1, 3.2 Myr), |δB| and |B R | are maximized. We show the quasi-periodic spatial changing by the online movies of Figure 4. Steady-state behavior is also seen from examining the mass flux inside 0.9 pc for t > 1.75 Myr. The radial profile of the net mass accretion rate appears to be roughly constant in time, which has been observed in simulation studies (e.g. Stone & Pringle 2001;Jiang et al. 2019). We took the time averages over 1.6 Myr, which is the longest time scale in the periodic oscillation of |B ϕ |. The time variation of the mass accretion rate due to MRI turbulence has been discussed in O'Neill et al. (2011), Hawley et al. (2011 and Hogg & Reynolds (2016). Since the mass flux varies with the averaging time interval, the time domain 0.744 < t < 1.75 Myr is not steady.

Direction Reversal and Vertical Transport
The periodical cycle of the mean toroidal field energy is observed to have the pattern of a quasi-periodic direction reversal. Figure 6 (top), the so-called the butterfly diagram, is the space (z-direction)-time evolution of the mean toroidal field direction at R = 0.715 pc. It shows that the direction periodically changes at a given z, and it also implies that the B-field escapes from the mid-plane. Figure 6 (bottom) is for the mean plasma-β, β ≡ 2P g /|B| 2 . The disk surface traces the bounding surface for β ≤ 10. This implies that the magnetic field becomes amplified near the mid-plane and is transferred to the high latitudes of the disk. The quasi-steady radial mass flow indicates the saturation of the MRI-driven turbulence.
In the MHD+Cooling/Heating state (t > 0.744Myr), the lowest plasma-β is in the mid-plane, and mean field transport is slower than that of the adiabatic MHD. The direction of B ϕ in 1.20 < t Myr < 2.27 is vertically stratified on blue around the mid-plane and red above that. The vertical transport changes from slow to fast at time t = 1.74 Myr. Recall that the energy of B ϕ increases and decreases in Figure 5, and this time is the maximum point. The same phenomenon occurs on the next reversal where it is red around the mid-plane and blue above that. We have confirmed the cycle through the mean toroidal field direction reversal and escape from the disk.
The vertical escape of the mean field is described as the vertical magnetic energy transport using the Poynting flux F P, z , where a rise speed v B denotes the vertical derivative of the pattern of the mean toroidal field (Salvesen et al. 2016b). The large rise speeds occur at the direction reversal, and large fall speeds appear around the mid-plane and |z| > 0.5 pc (Figure 7). The maximum rise speed v B is 37.3 km s −1 in |z| < 0.8 pc. In the same way, Alfvén speed and vertical speed are measured as |v A | ∼ 41.6 km s −1 and |v z | ∼ 44.0 km s −1 , and these speeds are comparable. The rise speed cannot reach the gravitational escape speed, 244.8 km s −1 at R = 0.715 pc.

Turbulent Velocity Field and Magnetoconvective Instability
In this subsection, we investigate the physical origin of the turbulent magnetic field shown in Section 3.1 and how it is maintained. Figure 8 shows the turbulent velocity field of the gas at t = 2.353 Myr. The maximum up-flow (blue arrows) is about 18% of the escape velocity at that position. The gas circulates with downward and upward flows in the ambient, where log {P g /(γ g − 1)} −4. The velocity inside the disk is relatively smaller than that in the ambient. The magnitude of the turbulent velocity for R 0.9 pc is comparable to the sound speed. That for R 0.9 pc is roughly 10% of the sound speed. The direction of vectors does not coincide with the turbulent and mean magnetic fields; however, part of the gas falls into the mid-plane along the magnetic field lines. Each component of the turbulent field and velocity has a difference only within one order of magnitude. Figure 9 shows the vertical structures of the magnetic, thermal, and kinetic energies at a given radius at t = 2.533 Myr. Around the mid-plane, the mean toroidal magnetic field dominates the turbulent magnetic and kinetic energies (ρ|δv| 2 /2), both of which are comparable to the thermal energy (P g /(γ g − 1)). The mean field decreases more rapidly than the turbulent component with z. As a result, there are regions where δB 2 ϕ > B 2 ϕ (0.2 |z| pc 0.8). The transition between the turbulent component and the mean field also occurs at the disk surface (|z| ∼ 0.05 pc), where the mean field reverses its direction as seen in Figure 6. At higher latitude (|z| 0.3 pc), contrary to the mid-plane (z = 0 pc), turbulent toroidal field energy (δB 2 ϕ /2) is about 10 % of the mean thermal energy, but it is comparable to or a few times larger than the kinetic energy.
The vertical random motion could be related to the magnetoconvective instability (or interchange instability, see e.g., Acheson 1979). The unstable criterion is d dz This criterion corresponds to the convective instability for the gas with decreasing specific entropy s with large |z|, i.e.
Replacing P g and γ g in the relation (27) with P B and γ B , where the magnetic pressure P B = |B| 2 /2, the criterion (26) is obtained. Here γ B = 1 + P B /U B = 2 for the energy density U B = |B| 2 /2 (e.g., Kulsrud 2005). The criterion (26) can be written as d (|B ϕ |/ρ) /dz < 0. Therefore, the instability occurs when the mass frozen in the magnetic flux tube per unit length increases with increasing |z|. One should note that the criteria (27) and (26) are necessary conditions for convection. We compare the two criteria with the vertical variation of the modified plasmaβ, i.e.β ≡ P g + P K /P B in Figure 10, where P K ≡ ρδv 2 . Note that the strong ram pressure exists around the mid-plane, and P K ∼ P g , but magnetic energy dominates. In the bottom panel of Figure 10, hydrodynamic convection is essentially stable, i.e. the positive entropy gradient shown in shadows in red. In contrast, as seen in the top panel of Figure 10, there are many unstable regions for the magnetoconvective instability (blue shadows). From the conditions, i.e. d dz PB ρ 2 < 0(unstable), d dz Pg ρ γg > 0(stable), after some algebraic calculations to vanish the density gradient, we can derive When d|B ϕ |/dz ∼ 0 (i.e. when field direction reversal occurs) and/or γ g /2 ∼ 1, this condition satisfies dβ/dz > 0. Figure 10 shows that the regions with dβ/dz > 0 correspond to unstable regions for magnetoconvective instability. The two panels of Figure 11 show the unstable regions for the two criteria, (26) and (27). It shows that the unstable regions for the magnetoconvection form belt-like layers, where the mean toroidal field direction reverses (see, B ϕ on R-z plane in Figure 4). The unstable layers move with the rising direction pattern of the mean toroidal field. The width of the unstable layers remains approximately constant around 0.05-0.1 pc. The bottom panel implies that the system is convective stable. . Vertical distribution of mean energy densities in MHD+Cooling/Heating state at t = 2.353 Myr. The mean energy density in the turbulent fields δB 2 i and ρδv 2 i (i = R, ϕ, z) are averaged over the azimuthal direction. Normalized energy densities are the rotational energy, ρ0vK /2, estimated by substituting the number density n = 10 2 cm −3 and Keplerian rotation vK at R = 0.715 pc. Colors are the thermal energy of green, the turbulent kinetic energy of black, and the toroidal magnetic energy of blue, respectively. Dotted and solid lines denote the turbulent and mean component.

Thermal state on the magnetic activity
Criterion (27) Criterion (26) for gas pressure for magnetic pressure Figure 10. Vertical distribution of the mean modified plasma-β (solid line) and the criterion of convective instabilities (shadows) at t = 2.353 Myr. The mean modified plasma-β represents Pg + PK /PB, where PK is the ram pressure evaluated by the turbulent kinetic energy, ρ|δv| 2 . The vertical distribution stability is shaded in red, and the unstability in blue. Stability evaluation is performed using mean quantities gradient for Equation 26 (upper) and Equation 27 (lower). Figure 12 shows the vertical distributions of the density and temperature at R = 0.715 pc at t = 2.353 Myr. For |z| < 0.05 pc, the gas is cold (∼ 8 × 10 2 K) and dense (n ∼ 10 3 cm −3 ). As seen in Figure 2(c), this cold disk extends to R ∼ 2 pc. Outside this cold, dense disk, the gas is warm (∼ 10 4 K) and less dense (n ∼ 10 cm −3 ); therefore, they are roughly in pressure equilibrium. The temperature of the ambient region (|z| > 0.25 pc) increases continuously until ∼ 10 6 K. Figure 12 also shows that there is a large azimuthal fluctuation around the mean values. The density fluctuation is in the order of unity or less while the temperature, as a maximum, fluctuates by three orders of magnitude. For the cold disk (|z| < 0.05 pc), the minimum temperature reaches the lower limit, and the maximum does not exceed 10 4 K of the warm disk.
To quantify the temperature and density fluctuations in phase space, we plot the gas mass fraction as a function of gas pressure and number density in the top panel of Figure 13. A large amount of gas mass is collected in characteristic regions over the gas mass fraction log(M/M tot ) > −5. These are multi-phase states created by the thermal instability. The gas is also distributed in a wide range of densities and temperatures, but a large fraction of the gas is in a state with P g /k B 10 5 and n 1. The bottom panel in Figure 13 shows the spatial distributions of four thermal states shown in the top panel on the R-z plane. The colors represent the region enclosed by the same colors on P g -n plane.

Unstable Stable
Unstable Stable The low-temperature gases is patchy in the warm (∼ 10 4 K) thick disk shown in blue. The mass fraction of this low temperature gas is small. Although the mass fraction of gas in the red region is small, The green region refers to the dense gas in a cold (∼ 100 K) phase. The gas shown in yellow in the top panel mostly forms the thick disk at R > 0.9 pc as shown in the bottom panel, where they are approximately isothermal (therefore, the gas pressure is constant). Figure 13. Top panel is a 2D histogram of the gas pressure in the vertical axis and the number density in the vertical axis at t = 2.352 Myr. Color contour denotes the mass fraction distribution occupying cells of ∆(log n) = ∆(log Pg/kB) = 0.01 and averaged over the number of mesh points in each Pg-n space cell. Bottom panel is the R-z plane distribution corresponding to some regions enclosing each colors on the top panel; Light blue: 3.4 < log(T [K])< 4.6, −1 < log(n[cm −3 ]) < 2, Red: 2.0 < log T < 3.4, 3 < log(Pg/kB [Kcm −3 ]) < 5, Yellow: 1 < log T < 3.4, 5 < log(Pg/kB) < 5.8, Green: 1 < log T < 3.6, 5.8 < log(Pg/kB ) < 6.8. Note that the gas temperature denotes a constant slope as d log(Pg/kB)/d log n, assuming the ideal gas.

Gas mass fraction
In Figure 14, density and temperature for R < 0.9 pc (i.e., the quasi-steady cold, thin gas disk) are plotted as a function of the magnetic pressure P B and the thermal pressure P g at t = 2.35 Myr. In the region where magnetic pressure dominates (β < 1), the number density of the left panel is distributed in a wide range from n ∼ 10 −2 to 10 4 cm −3 for a given gas pressure. It shows that the temperature is lower than T ∼ 10 5 K, and the strongly magnetized gas consists of cold gas with T < 1000 K, which dominates the total mass ( Figure 13). In the high-β domain, the temperature and density fixed a gas pressure are not sensitive to changes in magnetic field. The thermal state is determined by the compression and expansion of the gas pressure.

Direction reversal in the low β MRI
We found that the low β disk formed by radiative cooling and heating is discernible by the direction reversal (Fig.  6) and vertical transfer of the magnetic field (Fig. 7). A long reversal period are observed in the strong toroidal field with low β compared to the turbulence with high β. This trend is also found in the local 3D simulations of an isothermal gas (Bai & Stone 2013;Salvesen et al. 2016b). The mean plasma-β, β ∼ 0.1 − 0.4 in the mid-plane is also similar to our result (see, Fig. 10; Table 2 of Salvesen et al. 2016b). Compared to an analytical model (Begelman et al. 2015), it was demonstrated that heating efficiency, defined as the ratio between the dissipation rate (e.g. the magnetic reconnection) and mean toroidal field production rate, decreases with β. Salvesen et al. (2016b) showed that assuming steady Poynting flux, the period of the direction reversal is proportional to the rotational period (T cycle = 2π/Ωξ −1 B ) and is determined by the phenomenological parameters; ξ B = νη B /(β − ν + 1), where v B = η B Ωz (see also Equation 25) and ν is the degree of turbulent heating.
The turbulent heating caused by the nonlinear MRI is weakened by the small plasma-β and/or thermal pressure. Figure 15 shows that the turbulent magnetic energy (δB ϕ ; blue dotted line) cannot exceed the thermal energy (∼ P g ; green solid lines), i.e. P g /δB 2 ϕ > 10. Therefore, the MRI in the mid-plane is essential to amplify δB ϕ , since P g cannot be converted to δB ϕ . We rewrite plasma-β as, where assuming β −1 = (δB 2 ϕ /2 + B ϕ 2 /2)/P g . The plasma-β is decided by the thermal (P g /δB 2 ϕ ) and mean magnetic (B ϕ 2 /δB 2 ϕ ) energies. Given the ratio P g /δB 2 ϕ , the mean field dominant (B ϕ 2 ≥ δB 2 ϕ ) imposes an upper limit on the plasma-β, hence β ≤ 0.5P g /δB 2 ϕ . The lower one is limited by the MRI unstable condition, β > c s /v K (Begelman & Pringle 2007). Actually, the mean plasma β > 0.1 around the mid-plane in Figure 10 is consistent with this condition from the Keplerian rotation v K at R = 0.715 pc and mean temperature in Figure 12. We obtain the lower limit plasma β, We have showed by the simulations directly that as a consequence of satisfying this condition, the low β accounts for the mean field dominated disk. The field direction reversal as seen in Figure 6 is caused by the changes in balance between the turbulent energy and the mean field energy. When the direction of the mean field, which is represented by two shaded colors, reverses (e.g., t = 1.20, 2.12, 3.48 Myr in Figure 15), the energy density B 2 ϕ /2 (blue solid line) temporarily drops and the turbulent magnetic field (blue dotted line) is amplified. While the turbulent kinetic energy (black dotted line) exceeds the thermal energy, the turbulent magnetic energy (blue dotted line) is smaller than the thermal energy of cold gas (T ∼ 1000 K), i.e. β at this moment becomes large. The direction reversal is caused by a combination of MRI and the Parker instability ( Figure 16 or see e.g., Machida et al. 2013). A magnetic field line with a small radial fluctuation (panel a) is stretched Figure 15. Time variation of the energy densities in the mid-plane at R = 0.715 pc. Each energy density is normalized by the rotational energy with a density of ρ0 and Keplerian rotation at a radius of R = 0.715 pc. Solid and dotted lines represent the energy density of the mean and turbulent components, respectively. Different colors of solid lines denote the mean toroidal B-field energy (blue), the mean thermal energy (green), and the turbulent kinetic energy (gray). Shaded colors of red and blue signify the direction of the mean toroidal field. Figure 16. Schematic illustration of the MRI-Parker dynamo in the differential rotational disk. A magnetic field line is drawn as black curves with colored arrows. Red and blue arrows signify parallel and anti-parallel and the purple arrows indicate the rotational direction. (a) A magnetic field line is threaded by the toroidal field with a small radial turbulence on the R-ϕ plane at a certain height; (b) MRI and disk rotation stretch the magnetic field line to the radial and toroidal direction in R-ϕ plane; (c) Parker instability buoyantly escapes from the R-ϕ plane of (a).
to the radial and azimuthal directions (panel b). The stretched field line with a longer wavelength is selectively buoyed up toward higher latitudes by the Parker instability (panel c). As a result, for the R-ϕ plane, the buoyant field line has an opposite direction to the field line in the mid-plane. In the low β disk, direction reversal is possible by driving turbulence from the mean field; thus, the process to high β will be important (i.e. magnetoconvection in Section 3.3). Direction reversal was widely observed in adiabatic MHD simulations (e.g., O'Neill et al. 2011;Beckwith et al. 2011;Flock et al. 2012;Parkin & Bicknell 2013;Machida et al. 2013;Hogg & Reynolds 2016), but we found here for the first time that direction reversal also occurs in the low β disk with the multi-phase gas (10 < T [K]< 10 5 ).
It has been known that the MRI driven turbulence depends on numerical spatial resolution. Hawley et al. (2011Hawley et al. ( , 2013 applied the quality factor Q ϕ ≡ λ/(R∆ϕ), which is the resolution of the characteristic MRI wavelength (λ = 2πB ϕ / ρΩ 2 ). They identified the empirical condition as Q ϕ 20. However, for the toroidal field, the numerical convergence of adiabatic turbulence has not been clarified yet. The low β MRI in cold gas requires high resolution, where N ϕ = 2π/∆ϕ is the azimuthal resolution. In our model, N ϕ = 512 and we found that Q ϕ 20 (see, e.g. Kudoh & Wada 2018). Additionally, a comparison with N ϕ = 128 showed that the low resolution is Q ϕ < 10, and thus, there is no direction reversal in the cold and low β disk. Our simulations are in agreement with the estimation of Equation 31 and the empirical condition. In the long-term calculation, the turbulent magnetic field is sensitive to the numerical flux solver and the high-order accuracy. In order to reduce numerical dissipation, HLLD flux solver is employed (see, e.g. Hawley et al. 2013). Matsumoto et al. (2019, e.g., Section 3.2.1;4.3) pointed out that the high-order accuracy prevents the dissipation rather than the low-order scheme. We took a highly precise numerical approach, hence the differences from previous global 3D simulations.

Radiation pressure
AGNs emit enormous energy fluxes over a wide wavelength, and its feedback is important for the dynamics of the circumnuclear gas (e.g., Chan & Krolik 2016, 2017. However, IR radiative pressure is not dynamically effective (Namekata & Umemura 2016), and can not contribute to MRI in the cold gas. Notably, the anisotropic radiation pressure on the dust in the gas with the X-ray heating produces flows in a fountain-like manner in the central tens pc regions around the AGNs (Wada 2012). The radiation pressure is expected to be stronger than the magnetic pressure, i.e.
where κ d = 10 5 is the dust opacity and α d = 0.01 is the gas to dust ratio (Wada 2012). In a subsequent paper, we will investigate the effect of radiation pressure, how magnetic structures are changed, and how magnetic buoyancy or magnetoconvection can help.

SUMMARY
We studied the evolution of a magnetized multi-phase gas using global 3D MHD simulations in the pc-scale galactic nuclei. The simulation starts from an adiabatic state (L = 0 in Equation 5) with a weak toroidal field, β = 100, until the MHD turbulence is fully developed (β ∼ 0.6) for ∼ 25 rotational periods at R = 1 pc. Thereafter, the effects of the radiative cooling and the X-ray heating from the accretion disk around the SMBH are taken into account for additional 97 rotational periods (2.89 Myr). The magnetic pressure dominated disk is formed due to MRI cooperating with the radiative cooling. The quasi-steady state in a time sufficiently longer than the dynamical timescale is attained in the radii R < 0.9 pc for t > 1.75 Myr, as confirmed by a time variation of the mean toroidal field and a constant of the accretion rate averaged over variation timescale.
Major findings are: (i) The cold (< 10 3 K) gas forms a geometrically thin disk around the mid-plane. The warm (∼ 10 4 K) gas forms a thicker disk, and the hot (∼ 10 6 K) gas is distributed to higher latitudes (Figures 12, 13). The mass fraction (> 10 −5 ) on the P g -n plane is in the warm and cold phases ( Figure 13). The magnetic pressure is stronger (i.e. β(< 1)) in the cold, dense gas ( Figure 14).
(ii) The mean magnetic field is dominated by a toroidal component, and a strongly magnetized cold disk with β ∼ 0.02 is formed (Figures 2, 3, and 4). The mean toroidal field moves with the cold gas to a radial and vertical direction. The energy of the turbulent field is suppressed by the cooling effect; however, it is always comparable to or smaller than the thermal energy ( Figures 5, 15).
(iii) The turbulent motion in the multi-phase gas is observed in the R-z plane (Figure 8). Magnetoconvective instability plays a key role in maintaining turbulence for a long period,. The unstable condition d(|B ϕ |/ρ)/dz < 0 (Equation 26) coincides with the region in the modified plasma beta increasing vertically upward ( Figure 9) and in the belt-like mean field reversal in the R-z plane (Figure 11). The transition between the turbulent component and the mean field also occurs at the disk surface (|z| ∼ 0.05 pc and 0.2 |z| pc 0.8), where the mean field spatially reverses its direction (Figures 4, 9).
(iv) The quasi-steady state differs for the plasma β inside the disk. The high β disk is achieved by the saturation of |δB| 2 and the oscillation of B 2 ϕ ( Figure 5). In our simulations starting from the initial weak toroidal field (β = 100), the magnetic field strength amplified by the MRI remains β 1 and the oscillation timescale is about 10 rotational periods, as mentioned in the previous studies. The quasi-steady state in the low β disk is obtained by the periodic change of B 2 ϕ . The period in the low β state is more than 5 times longer (about 50 rotational periods) than that found in the high β state.
(v) We found that even for low β (∼ 0.02), the mean toroidal field shows direction reversals with time ( Figures 5,  6). This is caused by the transportation of the magnetic field vertically due to Parker instability (Figure 16), similar to the adiabatic state. The direction reversal of the mean B ϕ occurs, when the turbulent magnetic energy becomes larger than mean magnetic energy ( Figure 15). The moving speed of the magnetic field (Equation 25) estimated by the Poynting vertical flux is about 10% of the rotation speed. This speed becomes a maximum where the mean toroidal field direction reverses (Figure 7). To continue this cycle, mean magnetic flux transport from mid-plane to vertical direction is important.

ADDITIONAL LINK
Movies of snapshots of Figure 2 and 4 are available in the following link https://astrophysics.jp/MHD torus/.

ACKNOWLEDGMENTS
The authors are grateful to the anonymous referee for constructive comments and suggestions. We thank Ryoji Matsumoto, Mami Machida, and Yusuke Tsukamoto for their constructive comments and discussion, and the CANS+ developer team for the numerical techniques. Numerical computations were performed on Cray XC50 and XC30 systems at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI Grant Number 16H03959.

A. NON-IDEAL MHD EFFECTS
In this paper, we solved ideal MHD equations, whereby the gas in the circumnuclear region moves together with the magnetic field. This would be justified because we assumed that the ionization degree is (x ≡ n i /n n ), x ∼ 10 −4 following Meijerink & Spaans (2005). However, this assumption is incorrect if the magnetic Reynolds numbers for the dissipation processes in terms of the ohmic effect, Hall effect, and ambipolar diffusion are smaller than unity. Here, we confirm this.
The electric field in the induction equation (Equation 2) is replaced by the generalized Ohm's law (e.g., Braginskii 1965), where J ⊥ is the current perpendicular to the magnetic field. The three terms in the r.h.s of Equation (A1) are the magnetic dissipation of the ohmic, Hall, and ambipolar terms, respectively. The coefficients are formulated as follows: η = m e (ν ei + ν en ) 4πe 2 n e , η A = ρ n ρ 2 |B| 2 4π (ρ i ν in + ρ e ν en ) , where indices i, e, n denote the particle species of ion, electron, and neutral hydrogen, respectively. Here, ν ab denotes the collisional frequency of a particle "a" with a particle "b". Collisional frequencies are given by Spitzer (1962) assuming elastic collision, where the Coulomb logarithm ln Λ is about one order of magnitude. The electron-neutron collision is not effective, ν en /ν ei ∼ 0.52/ ln Λ(10 4 /x)(T /10 3 K) 2 < 1. In dense gas, ν in = 3.5 × 10 13 ρ n is often used (Draine et al. 1983). The magnetic Reynolds number (R M ≡ LV /η) is defined as the ratio of the |v × B| term to the dissipation term in Equation A1. We adopt the typical advection scale, L ∼ 1 pc and v ∼ v A ∼ 10 km s −1 , and the magnetic Reynolds numbers in each dissipation are, Equations (A8) and (A9) imply that we can ignore the ohmic and Hall dissipation. The ambipolar diffusion may be important for very strong magnetic fields (e.g. ≫ 1 mG) and/or diffuse media (n n < 100 cm −3 ), excepted for high ionization, e.g. x 10 −4 .