DYNAMO ACTIVITIES DRIVEN BY MAGNETOROTATIONAL INSTABILITY AND THE PARKER INSTABILITY IN GALACTIC GASEOUS DISKS

We solved the resistive MHD equations by using a modified Lax–Wendroff scheme with artificial viscosity (Rubin & Burstein 1967; Richtmyer & Morton 1967):

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot { \rho {\boldsymbol {v}}} = 0 \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial {\boldsymbol {v}}}{\partial t} + ({\boldsymbol {v}} \cdot \nabla){\boldsymbol {v}} = -\frac{1}{\rho } \nabla P - \nabla \phi + \frac{1}{4 \pi \rho } (\nabla \times {\boldsymbol {B}}) \times {\boldsymbol {B}} \end{equation} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\partial }{\partial t} \left(\frac{1}{2}\rho v^2 + \frac{B^2}{8 \pi } + \frac{P}{\gamma -1} \right)\nonumber\\ &&\qquad+ \nabla \cdot \left[ \left(\frac{1}{2}\rho {\boldsymbol {v}}^2 + \frac{\gamma P}{\gamma -1} \right) {\boldsymbol {v}} + \frac{c}{4 \pi } {\boldsymbol {E}} \times {\boldsymbol {B}} \right]= - \rho {\boldsymbol {v}} \nabla \phi \nonumber\\ \end{eqnarray} \tag{ 3 }$$

$$\begin{equation} \frac{\partial {\boldsymbol {B}}}{\partial t} = \nabla \times ({\boldsymbol {v}} \times {\boldsymbol {B}} - \eta \nabla \times {\boldsymbol {B}}), \end{equation} \tag{ 4 }$$

where ρ, P, ${\boldsymbol {v}}$, ${\boldsymbol {B}}$, and γ are the density, gas pressure, velocity, magnetic field, and the specific heat ratio, respectively, and ϕ is the gravitational potential proposed by Miyamoto & Nagai (1975). The electric field ${\boldsymbol {E}}$ is related to the magnetic field ${\boldsymbol {B}}$ by Ohm's law. We assume the anomalous resistivity η

$$\begin{equation} \eta = \left\lbrace \begin{array}{@{}ll@{}} \eta _0 (v_{\rm d}/v_{\rm c} - 1)^2, & v_{\rm d} \ge v_{\rm c}, \\ 0,& v_{\rm d} \le v_{\rm c}, \end{array} \right. \end{equation} \tag{ 5 }$$

which sets in when the electron-ion drift speed v_d = j/ρ, where j is the current density, exceeds the critical velocity v_c (Yokoyama & Shibata 1994).

We assumed an equilibrium torus threaded by weak toroidal magnetic fields (Okada et al. 1989). The torus is embedded in a hot, non-rotating spherical corona. We assume that the torus has constant specific angular momentum L at ϖ = ϖ_b and assume a polytropic equation of state P = Kρ^γ, where K is constant. The density distribution is same as that for Equation (8) in Nishikori et al. (2006).

We consider a gas disk composed of one-component interstellar gas. Radiative cooling and self-gravity of gas are ignored in this paper.

We adopt a cylindrical coordinate system (ϖ, φ, z). The units of length and velocity in this paper are $\varpi _0 =1 \hspace{2.84526pt} \rm {kpc}$ and $v_0=(GM/\varpi _0)^{1/2} = 207 \hspace{2.84526pt} {{\rm km\,s^{-1}}}$, respectively, where M = 10¹⁰ M_☉. The unit time is t₀ = ϖ₀/v₀. Other numerical units are listed in Table 1.

The numbers of grid points we used are (N_ϖ, N_φ, N_z) = (250, 128, 640). The grid size is Δϖ/ϖ₀ = 0.05, Δz/ϖ₀ = 0.01 for 0  ⩽  ϖ/ϖ₀ ⩽ 6.0, and |z|/ϖ₀ ⩽ 2.0, and otherwise they increase with ϖ and z, respectively. The grid size in the azimuthal direction is Δφ = 2π/128. The outer boundaries at ϖ = 56ϖ₀ and |z| = 10ϖ₀ are free boundaries where waves can be transmitted. We applied periodic boundary conditions in the azimuthal direction. An inner absorbing boundary condition is imposed at $r = \sqrt{\varpi ^2 + z^2} = r_{\rm in} = 0.8 \varpi _0$, since the gas accreted onto the central region of the galaxy will be converted to stars or swallowed by the black hole. To illuminate the effect of the boundary condition at the equatorial plane, we calculate two models: one is model PSYM, which imposes a symmetric boundary condition at the equatorial plane. This model is identical to that in Nishikori et al. (2006). The other is model ASYM, in which we include the equatorial plane inside the simulation region, and imposed no boundary condition at the equatorial plane. In this paper, we adopted model parameters ϖ_b = 10ϖ₀, β_b = 100 at (ϖ, z) = (ϖ_b, 0), γ = 5/3, v_c = 100v₀, and η₀ = 0.1v₀ϖ₀. The thermal energy of the torus is parameterized by E_th = c²_sb/(γv²₀) = Kρ_b^{γ − 1}/v²₀ = 0.05 where c_sb and ρ_b are the sound speed and density at (ϖ, z) = (ϖ_b, 0), respectively. The initial temperature of the torus is T ∼ 2 × 10⁵ K. This mildly hot plasma mimics the mixture of hot (T ∼ 10⁶ K) and warm (T ∼ 10⁴ K) components of interstellar matter, which occupy a large volume of the galactic disks.

Figure 1(a) shows the time evolution of the mass accretion rate measured at ϖ/ϖ₀ = 2.5. After t/t₀ ∼ 200, the gas starts to inflow and the mass accretion rate increases linearly. Subsequently (t/t₀ > 600) the accretion rate saturates and becomes roughly constant. The evolution of the mass accretion rate in model ASYM (black) is similar to that in model PSYM (gray), where the mass accretion rate for model PSYM is doubled since model PSYM includes only the region above the equatorial plane.

Zoom In Zoom Out Reset image size

Figures 1(b)–(d) show the time evolution of the plasma β and the magnetic energy of each component averaged in the region where 2 ⩽ ϖ/ϖ₀ ⩽ 5, |z|/ϖ₀ ⩽ 1, and 0 ⩽ φ ⩽ 2π, respectively. It is clear that plasma β decreases in the linear stage as the magnetic energy increases. In both models, the averaged magnetic energy first increases exponentially. After that, it saturates and becomes roughly constant. The time evolution of magnetic fields is similar between model ASYM and model PSYM. The time evolutions of the growth and saturation in the mass accretion rate are quite similar to those of the averaged magnetic energy. This means that the magnetic turbulence driven by the MRI is the main cause of the angular momentum transport which drives mass accretion. In the remaining part of this section, we discuss the results of model ASYM.

In order to check the magnetic field structure, we analyze the time evolution of mean azimuthal magnetic fields. The mean fields are computed by the same method as Nishikori et al. (2006). Figure 2(a) shows the time evolution of the mean azimuthal magnetic fields. The black curve shows B_φ averaged in the region where 5 < ϖ/ϖ₀ < 6, 0 ⩽ φ ⩽ 2π, and 0 < z/ϖ₀ < 1, and the gray curve shows B_φ averaged in the region where 5 < ϖ/ϖ₀ < 6, 0 ⩽ φ ⩽ 2π, and 1 < z/ϖ₀ < 3. Black and gray curves correspond to the disk region and halo region, respectively. The azimuthal magnetic fields reverse their direction on a timescale t ∼ 300t₀ ∼ 1.5 Gyr. This timescale is comparable to that of the buoyant rise of azimuthal magnetic fields. The halo fields also reverse their directions on the same timescale as the disk component.

Zoom In Zoom Out Reset image size

The radial distribution of the mean azimuthal magnetic fields at t = 1000 t₀ are shown in Figure 2(b). Black and gray curves show the averages in the disk (0 ⩽ φ ⩽ 2π and 0 < z/ϖ₀ < 1) and in the halo (0 ⩽ φ ⩽ 2π and 1 < z/ϖ₀ < 3), respectively. Since the magnetic fields in the disk become turbulent, the direction of azimuthal magnetic fields changes frequently near the equatorial plane. On the other hand, long wavelength magnetic loops are formed in the halo region.

Figure 3(a) shows the radial distribution of the specific angular momentum at the equatorial plane. Although the rotation is assumed only inside the initial torus (gray curve in Figure 3(a)), the torus deforms its shape from a torus to a disk by redistributing angular momentum (black curve in Figure 3(a)). Since the outflow created by the MHD-driven dynamo transports angular momentum toward the halo, the halo begins to rotate differentially, as can be seen in isocontours of the azimuthally averaged rotation speed (Figure 3(b)). Figures 3(a) and (b) show that the azimuthal velocity v_φ is almost constant and close to v₀ in the equatorial plane. The rotation period at ϖ/ϖ₀ = a is P_rot ∼ 2πa(ϖ₀/v₀) = 2πat₀.

Zoom In Zoom Out Reset image size

We show the 3D structure of the mean magnetic fields in Figure 4. The colored surfaces represent isodensity surfaces, and the curves show magnetic field lines. Light brown curves show field lines passing through the plane at z/ϖ₀ = 5 and ϖ/ϖ₀ < 1. The vertical magnetic fields around the rotation axis are produced by the magnetic pressure-driven outflow near the galactic center. Color on the curves depicts the direction of the azimuthal magnetic fields. Red curves show a positive direction of the azimuthal magnetic fields (counterclockwise direction) and blue shows a negative direction (clockwise). The figure indicates that the azimuthal field reverses the direction with the radius.

Zoom In Zoom Out Reset image size

Figure 4 indicates that long wavelength magnetic loops are formed around the disk–halo interface. Therefore, we discuss the formation mechanism of the magnetic loops. Horizontal magnetic fields embedded in a gravitationally stratified atmosphere become unstable against long wavelength undular perturbations. This undular mode of magnetic buoyancy instability, called the Parker instability, creates buoyantly rising magnetic loops (Parker 1966). Nonlinear growth of the Parker instability in gravitationally stratified disks was studied by MHD simulations by Matsumoto et al. (1988). They showed that magnetic loops continue to rise when β < 5. In a weakly magnetized region where β > 5, the Parker instability only drives nonlinear oscillations (Matsumoto et al. 1990).

In differentially rotating disks, the MRI couples with the Parker instability. The growth timescale of the non-axisymmetric the MRI is t_MRI ∼ 1/(0.1Ω) ∼ 10H/c_s where c_s is the sound speed. On the other hand, the growth timescale of the Parker instability is $t_{\rm PI} \sim 5H/v_{\rm A} \sim 5\sqrt{\beta \gamma /2} (H/c_{\rm s})$ where v_A and γ are the Alfvén speed and the specific heat ratio, respectively. This timescale becomes comparable to the growth time of the MRI when β ∼ 5. Since magnetic turbulence produced by the MRI limits the horizontal length of coherent magnetic fields, the growth of the Parker instability can be suppressed in a weakly magnetized region inside the disk. On the other hand, in regions where β < 5, the nonlinear growth of the Parker instability creates magnetic loops buoyantly rising to the disk corona (Machida et al. 2000, 2009).

Figure 5(a) shows an example of a magnetic loop at t = 1000 t₀. Magnetic loops are identified by the same algorithm as that reported in Machida et al. (2009). Figures 5(b)–(d) show the distribution of the vertical velocity, plasma β, and density, respectively, along the magnetic field line depicted in Figure 5(a). A black circle shows the starting point of the integration of a magnetic field line. The positive vertical velocity indicates that the magnetic loop is rising. The value of plasma β decreases from the footpoint of the loop where β ∼ 100 to the loop top where β < 5 (Figure 5(d)). The density also decreases from the footpoints to the loop top. These are typical structures of buoyantly rising magnetic loops formed by the Parker instability.

Zoom In Zoom Out Reset image size

From the numerical results presented so far, the initial magnetic field was found to be amplified and transported to the halo region after several rotation periods. In order to investigate the time evolution of the azimuthal magnetic field, we show a butterfly diagram of the azimuthal magnetic fields at ϖ/ϖ₀ = 2 (Figure 6(a)). The horizontal axis shows the time t/t₀, and the vertical axis shows the height z/ϖ₀. Color denotes the direction of the azimuthal magnetic fields. The white curve shows the isocontour where β = 5 above the equatorial plane. We should note that there is no magnetic flux at ϖ/ϖ₀ = 2 initially. As the time elapses, the accretion of the interstellar gas driven by the MRI transports magnetic fields from ϖ/ϖ₀ = 10 to ϖ/ϖ₀ = 2. The amplified azimuthal magnetic fields change their directions quasi-periodically. Turbulent magnetic fields are dominant around the equatorial plane.

Zoom In Zoom Out Reset image size

When the magnetic tension becomes comparable to that of the turbulent motion in the disk surface and the coherent length of the magnetic field line increases, the magnetic flux buoyantly escapes from the disk due to the Parker instability. We find that the instability significantly grows, when the plasma β decreases to around 5. The averaged escape velocity of the flux is about 25 km s⁻¹, approximately equal to the Alfvén velocity.

The periodic reversals of azimuthal magnetic fields can also be seen in the correlation between the fields above and below the equatorial plane (Figure 6(b)). White shows the region where B_φ(z) · B_φ(− z) > 0, and black shows the region where B_φ(z) · B_φ(− z) < 0. We can see fine structures of turbulence near the equatorial plane. In the halo region where z/ϖ₀ > 1, the color of the outgoing flux changes quasi-periodically. This means that the disk magnetic topology changes between the symmetric state and the anti-symmetric state.

Finally, we check the spatial distribution of the azimuthal magnetic fields. The azimuthal magnetic fields in the ϖ − z plane averaged in the region 0 ⩽ φ ⩽ 2π at t/t₀ = 1000 is shown in Figure 7(a). We can see that a few μG fields are distributed as high as at z/ϖ₀ = 5. Figures 7(b) and (c) show the distribution of the azimuthal magnetic fields in the ϖ − φ plane, where the azimuthal magnetic fields are averaged in the halo (1 < z/ϖ₀ < 1.5) and in the disk (0 < z/ϖ₀ < 0.5), respectively. Since turbulent magnetic fields are dominant inside the disk, short magnetic filaments are formed (Figure 7(c)). Long and strong azimuthal magnetic sectors are formed in the halo region because the plasma β in the halo region becomes lower than that in the disk due to the density decrease (Figure 7(b)). Therefore, the magnetic tension of buoyantly rising magnetic loop exceeds the ram pressure by the turbulent motion around the disk surface and in the halo. These magnetic loops buoyantly escape from the disk by the Parker instability and create the rising magnetic fluxes in the butterfly diagram shown in Figure 6(a).

Zoom In Zoom Out Reset image size

In order to study the dependence of the numerical results on azimuthal resolution, we carried out simulations in which we used 64 (model ASYM64) and 256 (model ASYM256) grid points in 0 ⩽ φ ⩽ 2π. Other parameters are the same as those in model ASYM.

Figure 8(a) shows the time evolution of the magnetic energy averaged in the same region as Figure 1(c). Black, dashed, and gray curves show models ASYM, ASYM256, and ASYM64, respectively. This panel indicates that the saturation level of the magnetic energy near the equatorial plane slightly decreases as the resolution increases because shorter wavelength turbulence inside the disk can be resolved better in the model with higher azimuthal resolution, so that more magnetic energy dissipates in the nonlinear stage. The time evolution of the azimuthal magnetic fields averaged in the region where 5 < ϖ/ϖ₀ < 6, 0 < z/ϖ₀ < 1, and 0 ⩽ φ ⩽ 2π is shown in Figure 8(b). Colors are the same as in Figure 8(a). The amplitude of the azimuthal magnetic fields and the timescale of the field reversal are almost similar in models ASYM and ASYM256. In order to resolve the most unstable wavelength of the MRI, we need more than 20 grid points per disk thickness for azimuthal direction (Hawley et al. 2011). Although ASYM does not satisfy this condition, the qualitative behavior of the results for ASYM is identical to that of ASYM256, which satisfies the condition (see Figure 8(b)). For example, the ratio of the radial magnetic energy to the azimuthal magnetic energy (B²_ϖ/B_φ²) is 0.12 in model ASYM and 0.13 in model ASYM256. Hence, we conclude that the numerical results of model PSYM and ASYM in this paper obtained by simulations with 128 grid points in the azimuthal direction do not significantly differ from results with higher azimuthal resolution in the nonlinear stage.

Zoom In Zoom Out Reset image size