# GLOBAL SIMULATIONS OF DIFFERENTIALLY ROTATING MAGNETIZED DISKS: FORMATION OF LOW-BETA FILAMENTS AND STRUCTURED CORONA

## 1 INTRODUCTION

Magnetic fields in differentially rotating disks play essential roles in the angular momentum transport which enable the accretion and various activities such as X-ray flares and jet formation. Motivated by the Skylab observations of the solar corona, Galeev, Rosner, & Vaiana (1979) proposed a model of magnetically structured corona in accretion disks consisting of X-ray emitting magnetic loops. The magnetic loops can be created due to the buoyant rise of magnetic flux tubes (or flux sheet) from the interior of the accretion disk. Matsumoto et al. (1988) carried out two-dimensional magnetohydrodynamic (MHD) simulations of the Parker instability (Parker 1966) in nonuniform gravitational fields which mimic those in accretion disks. They showed that when $`\beta =P_{\mathrm{gas}}/P_{\mathrm{mag}}1`$, a plane-parallel disk deforms itself into evacuated undulating magnetic loops and dense blobs accumulated in the valley of magnetic field lines. The effects of rotation and shear flow, however, were not included in their simulations.

Shibata, Tajima, & Matsumoto (1990) carried out two-dimensional MHD simulations of the Parker instability including the effects of shear flow and suggested that magnetic accretion disks fall into two types; gas pressure dominated (high-$`\beta `$) solar type disks, and magnetic pressure dominated (low-$`\beta `$) cataclysmic disks. In high-$`\beta `$ ($`\beta 1`$) disks, magnetic flux escapes from the disk more efficiently as $`\beta `$ decreases.

Several authors (Hawley, Gammie & Balbus 1995 ; Matsumoto & Tajima 1995 ; Brandenburg et al. 1995; Stone et al. 1996) have reported the results of three-dimensional local MHD simulations of magnetized accretion disks by adopting a shearing box approximation (Hawley et al. 1995). In differentially rotating disks, magnetorotational instability (Balbus & Hawley 1991) couples with the Parker instability (Foglizzo & Tagger 1995). Since Parker instability grows for long wave length perturbations along magnetic field lines, non-local effects may affect the stability and nonlinear evolution.

Results of global 3D MHD simulations including vertical gravity were reported by Matsumoto (1999) and Hawley (1999). By adopting an initially gas pressure dominated torus threaded by toroidal magnetic fields, Matsumoto (1999) showed that magnetic energy is amplified exponentially owing to the growth of the Balbus & Hawley instability and that the system approaches a quasi-steady state with $`\beta 10`$. Matsumoto (1999) assumed $`\beta =constant`$ in the torus at the initial state. When $`\beta _01`$, the deviation from magneto-rotational equilibrium introduces large amplitude perturbations.

In this paper, we present the results of 3D MHD simulations starting from an equilibrium MHD torus threaded by initially equipartition strength ($`\beta 1`$) toroidal magnetic fields.

## 2 MODELS AND NUMERICAL METHODS

The basic equations we use are ideal MHD equations in cylindrical coordinate system $`(r,\varphi ,z)`$. We assume that the gas is inviscid and evolves adiabatically. Since we neglect radiative cooling, our numerical simulations postulates that the disk is advection-dominated (see Kato, Fukue, & Mineshige 1998 for a review).

The initial condition is an equilibrium model of an axisymmetric MHD torus threaded by toroidal magnetic fields (Okada et al. 1989). We assume that the torus is embedded in hot, isothermal, non-rotating, spherical coronal gas. For gravity, we use the Newtonian potential. We neglect the self-gravity of the gas. At the initial state the torus is assumed to have a constant specific angular momentum $`L`$.

We assume polytropic equation of state $`P=K\rho ^\gamma `$ where $`K`$ is a constant, and $`\gamma `$ is the specific heat ratio. According to Okada et al. (1989), we assume

$$v_A^2=\frac{B_\varphi ^2}{4\pi \rho }=H(\rho r^2)^{\gamma 1}$$
(1)   
where $`v_A`$ is the Alfvén speed and $`H`$ is a constant. For normalization we take the radius $`r_0`$ where the rotation speed $`L/r_0`$ equals the Keplerian velocity $`v_{K0}=(GM/r_0)^{1/2}`$ as unit radius and set $`\rho _0=v_{K0}=1`$ at $`r=r_0`$. Using these normalizations, we can integrate the equation of motion into the potential form;

$$\mathrm{\Psi }=\frac{1}{R}+\frac{L^2}{2r^2}+\frac{1}{\gamma 1}v_s^2+\frac{\gamma }{2(\gamma 1)}v_A^2=\mathrm{\Psi }_0=constant,$$
(2)   
where $`v_s^2`$ is the square of the sound speed, $`R=(r^2+z^2)^{1/2}`$ and $`\mathrm{\Psi }_0=\mathrm{\Psi }(r_0,0)`$. By using equation (2), we obtain the density distribution.

$$\rho =\left\{\frac{\mathrm{max}[\mathrm{\Psi }_0+1/RL^2/(2r^2),0]}{K[\gamma /(\gamma 1)][1+\beta _0^1r^{2(\gamma 1)}]}\right\}^{1/(\gamma 1)}$$
(3)   
where $`\beta _02K/H`$ is the plasma $`\beta `$ at $`(r,z)=(r_0,0)`$. The parameters describing the structure of the MHD torus are $`\gamma `$, $`\beta _0`$, $`L`$ and $`K`$. In this paper we report the results of simulations for parameters $`\beta _0=1`$, $`\gamma =5/3`$, $`L=1`$, and $`K=0.05`$. The density of the halo at $`R=r_0`$ is taken to be $`\rho _{\mathrm{halo}}/\rho _0=10^3`$. The unit field strength is $`B_0=\rho _0^{1/2}v_{K0}`$.

We solve the ideal MHD equations in a cylindrical coordinate system by using a modified Lax-Wendroff scheme with artificial viscosity. We simulated only the upper half space ($`z0`$) and assumed that at the equatorial plane, $`\rho `$, $`v_r`$, $`v_\varphi `$, $`B_r`$, $`B_\varphi `$, and $`P`$ are symmetric and $`v_z`$ and $`B_z`$ are antisymmetric. The outer boundaries at $`r=r_{\mathrm{max}}`$ and at $`z=z_{\mathrm{max}}`$ are free boundaries where waves can transmit. The singularity at $`R=0`$ is treated by softening the gravitational potential near the gravitating center ($`R<0.2r_0`$). The number of grid points is $`(N_r,N_\varphi ,N_z)=(200,64,240)`$. To initiate non-axisymmetric evolution, small amplitude, random perturbations are imposed at $`t=0`$ for azimuthal velocity.

## 3 NUMERICAL RESULTS

Figure 1a shows the initial condition. Color scale denotes the density distribution and red curves depict magnetic field lines. Figure 1b shows density distribution and velocity vectors at $`t=6.2t_0`$ where $`t_0`$ is the orbit time $`t_0=2\pi r_0/v_{K0}`$. As the matter which lost angular momentum accretes to the central region, the MHD torus becomes disk-like. Velocity vectors indicate that matter flows out from the disk.

After the non-axisymmetric Balbus & Hawley instability grows, the inner region of the torus becomes turbulent. The turbulent motions tangle magnetic field lines in small scale and create numerous current sheets (or current filaments). Figure 1c shows the magnetic field lines projected onto the equatorial plane and the density distribution at $`z=0`$. Note that Figure 1c shows only the inner region where $`2.5x/r_02.5`$ and $`2.5y/r_02.5`$. In large scales, magnetic field lines and density distribution show low azimuthal wave number spiral structure. Figure 1d shows that magnetic loop structures similar to those in the solar corona are created. The yellow surfaces show strongly magnetized regions where $`|𝑩|=0.1B_0`$ and red curves show magnetic field lines at $`t=7.5t_0`$. The magnetic loops buoyantly rise from the disk due to the Parker instability. Numerical results indicate that small loops which newly appeared above the photosphere (emerging magnetic loops) develop into expanding coronal loops. We can observe magnetic loops elongated in the azimuthal direction and loops twisted by the rotation of the disk.

Figure 2 shows the isosurfaces of plasma-$`\beta `$ at $`t=6.2t_0`$. Orange surfaces show strongly magnetized regions where $`\beta =0.1`$. Yellow and green surfaces show the region where $`\beta =1`$ and $`\beta =10`$, respectively. Inside the disk, strongly magnetized, low-$`\beta `$ filaments are created. As we shall show below, low-$`\beta `$ region occupies a small fraction of the total volume. Intermittent magnetic structures (filamentary strongly magnetized regions) similar to those in the solar photosphere develop in the disk.

Figure 3a shows the spatial average of the magnetic energy, $`\mathrm{log}(B^2/8\pi /P)`$ (dashed curve), $`\mathrm{log}B^2/(8\pi P_0)`$ (dash-dotted curve) and $`\mathrm{log}(B_rB_\varphi /4\pi /P_0)`$ (solid curve) averaged in $`0.7r/r_01.3`$ and $`0z/r_00.3`$. This figure shows that magnetic energy in the equatorial region decreases within a few rotation period due to the buoyant escape of magnetic flux, and that averaged plasma $`\beta `$ oscillates quasi-periodically around $`\beta 5`$. The ratio of averaged Maxwell stress to initial equatorial pressure, $`\alpha _{th}`$ is $`\alpha _{th}0.07`$ when $`6t_0<t<12t_0`$. Figure 3b shows the time evolution of the accretion rate

$$\dot{M}_{\mathrm{acc}}=_0^{2\pi }_0^{0.3r_0}\rho v_rr𝑑z𝑑\phi $$
(4)   
at $`r=0.3r_0`$, and of the outflow rate

$$\dot{M}_{\mathrm{out}}=_0^{2\pi }_0^{r_{\mathrm{max}}}\rho v_zr𝑑r𝑑\phi $$
(5)   
at $`z=3.0r_0`$. Accretion rate increases and fluctuates around a mean value. Numerical results also indicate quasi-periodic ejection of the disk material. Figure 3c shows the Poynting flux $`(𝑬\times 𝑩/4\pi )_z`$ which passes through the plane at $`z=1.0,2.0,3.0r_0`$. After a few rotation period, magnetic flux is convected from the equatorial region to the disk surface. The volume filling factor of the region where $`\beta 0.3`$ is shown in Figure 3d. Strongly magnetized region where $`\beta 0.3`$ occupies about 8% of the total volume after a few rotation period. Later, the filling factor decreases to 2%. After about 10 rotation period, magnetic flux is regenerated in the disk and filling factor shows the second peak when the averaged $`1/\beta `$ is maximum. Following this second peak, Poynting flux at $`z=2r_0`$ increases, which suggest that magnetic flux escapes from the disk.

## 4 DISCUSSION

We showed by 3D global MHD simulations that when a differentially rotating torus is threaded by equipartition strength toroidal magnetic fields, magnetic loops emerging from the disk continue to rise and form coronal magnetic loops similar to those in the solar corona. Inside the disk, magnetic turbulence drives dynamo action which maintains magnetic fields and keeps the disk in a quasi-steady state with $`\beta 5`$. We successfully simulated the formation process of the solar type disk by 3D direct MHD simulations.

Numerical results indicate that in differentially rotating disks magnetic field lines globally show spiral structure with low azimuthal wave number but locally fluctuating components create numerous current sheets. We expect that magnetic reconnection taking place in current sheets generate hot plasmas which emit hard X-rays. When the disk is optically thin, such reconnection events may be observed as large amplitude sporadic X-ray time variations characteristic of low-states in black hole candidates (Kawaguchi et al. 1999).

We found that inside the torus, filamentary shaped, locally strongly magnetized, low-$`\beta `$ regions appear. Even when $`\beta 5`$ in average, low-$`\beta `$ regions where $`\beta 0.3`$ occupy 2-10% of the total volume. The low-$`\beta `$ filaments are embedded in high-$`\beta `$ plasma. Such an intermittent structure is common in magnetized astrophysical plasmas. Numerical results indicate that low-$`\beta `$ filaments are re-generated after they once disappear.

Although we assumed point gravity suitable for accretion disks, numerical results can qualitatively be applied to galactic gas disks. Our numerical results suggest that although $`\beta 1`$ in average, low-$`\beta `$ filaments exist in galactic gas disks. Magnetic reconnection taking place in such low-$`\beta `$ regions may heat the interstellar gas and create hot, X-ray emitting plasmas. Tanuma et al. (1999) proposed a model that magnetic reconnection in strongly magnetized ($`30\mu G`$) regions in our Galaxy creates hot plasma which emit 7KeV component of Galactic Ridge X-ray Emission (GRXE). The low-$`\beta `$ filaments can also confine the hot plasma and prevent it from escaping from the Galaxy.

We thank Drs. K. Shibata, T. Tajima, S. Mineshige, T. Kawaguchi and K. Makishima for discussion. Numerical computations were carried out by VPP300/16R at NAOJ. This work is supported in part by the Grant of Japan Science and Technology Corporation and the Grant-in-Aid of the Ministry of Education, Science, Sports, and Culture, Japan(07640348, 10147105).