Similarity between the Molecular Loops in the Galactic Center and the Solar Chromospheric Arch Filaments

We carried out two-dimensional magnetohydrodynamic simulations of the Galactic gas disk to show that the dense loop-like structures discovered by the Galactic center molecular cloud survey by NANTEN 4 m telescope can be formed by the buoyant rise of magnetic loops due to the Parker instability. At the initial state, we assumed a gravitationally stratified disk consisting of the cool layer ($T \sim 10^3$ K), warm layer ($T \sim 10^4$ K), and hot layer ($T \sim 10^5$ K). Simulation box is a local part of the disk containing the equatorial plane. The gravitational field is approximated by that of a point mass at the galactic center. The self-gravity, and the effects of the galactic rotation are ignored. Numerical results indicate that the length of the magnetic loops emerging from the disk is determined by the scale height of the hot layer ($\sim$ 100 pc at 1 kpc from the Galactic center). The loop length, velocity gradient along the loops and large velocity dispersions at their foot points are consistent with the NANTEN observations. We also show that the loops become top-heavy when the curvature of the loop is sufficiently small, so that the rising loop accumulates the overlying gas faster than sliding it down along the loop. This mechanism is similar to that of the formation of solar chromospheric arch filaments. The molecular loops emerge from the low temperature layer just like the dark filaments observed in the H$\alpha$ image of the emerging flux region of the sun.

(MHD) simulations of the Parker instability. They found that dense regions are formed in the valley of magnetic loops in the nonlinear stage of the Parker instability and that shock waves are formed at the footpoints of the rising magnetic loops where the gas infalling along the magnetic field lines collides with the dense gas near the equatorial plane of the disk. Shibata & Matsumoto (1991) applied the results of the MHD simulations to the formation of molecular clouds in the Orion region. Subsequently, Kamaya et al. (1996) studied the triggering of the Parker instability by supernova explosions. Basu, Mouschovias & Paleologou (1997) investigated the effect of the Parker instability on the structure of the interstellar medium, and Kim et al. (2000) applied it to galactic disks. Machida, Hayashi & Matsumoto (2000) reported the results of the three-dimensional (3D) global MHD simulations of the Parker instability in a differentially-rotating disk. They found that magnetic loops emerging from the disk form a structured corona above the disk. Fukui et al. (2006) presented the results of 2D MHD simulations of local volume of the Galactic disk. They suggested that the downflows can be the origin of the violent motion 1 and extensive heating of the molecular gas in the Galactic center. In their simulations, however, the cold component of the galactic gas that corresponds to molecular clouds was not taken into account and the mechanism of the formation of dense loop-like structures was not clear.
In the solar atmosphere, arch filaments are observed in the Hα images of the chromosphere. Isobe et al. (2006) found that dense filaments similar to Hα arch filaments are formed in the emerging flux. These filamentary structures are cool and dense above the chromosphere. This configuration in which the dense regions exist around the top of the emerging loops is a typical feature of the emerging loops and is a good approximation to what will be observed as dark features in the Hα image. They clarified the reason why the emerging loop becomes top-heavy on the basis of the results of 3D and 2D MHD simulations of the emergence of the magnetic loops from the convection zone, through the chromosphere, to the corona. The dense gas accumulated around the top of the rising loops fragments into filaments by the Rayleigh-Taylor instability, and slides down along the magnetic field lines. These filaments correspond to the arch filaments observed in the emerging flux regions of the sun.
In this paper, we present the results of 2D MHD simulations of the Galactic gas disk. We will consider a local volume of the Galactic gas disk that consists of a low temperature layer (the cold component of the Galactic gas disk), warm layer and hot layer through 2D MHD simulations. In section 2, we will present basic equations and the numerical model used in this study. The numerical results will be reported in section 3 and section 4 is devoted for discussion. Summary will be given in section 5.

Assumptions and Basic Equations
We adopt the local Cartesian coordinates (x, y, z) at radius r 0 from the Galactic center (see Figure 1). The x-direction is taken to be the azimuthal direction, and the z-direction is parallel to the rotational axis of the Galactic disk. The partial derivatives of the background medium with respect to the radial direction y are neglected (i.e., local approximation). The assumptions, basic equations and initial conditions are similar to those in Matsumoto et al. (1988). We will assume the following: (1) the medium is an ideal gas with the specific heat ratio γ = 1.05, (2) the viscosity and resistivity are neglected, (3) the effects of the rotation of the disk and self-gravity are neglected.
The basic equations are: and Here, g = (0, 0, g(z)) is the gravitational acceleration assumed to be where G is the gravitational constant. This vertical component of the gravitational acceleration is produced by the gravitational potential by a point mass and gives a simple but exact model for the accretion disk rotating around the point mass. In the case of galaxies, this is not an exact model, but at least the region of the equator can be approximated by this model, because the gravitational acceleration near the disk in galaxies is approximately proportional to the height from the equatorial plane of the disk. The radial component of gravity is assumed to be equal to the centrifugal force due to the rotation of the disk. The effects of the rotation will be discussed in section 4.4. Equations (1)-(4) are rendered dimension less by using normalizing constants r 0 , C S0 and ρ 0 , where C S0 is the sound speed in the mid-temperature region and ρ 0 is the unperturbed density at the equatorial plane. When the numerical results are compared with observations, we use the units of length, velocity, and time in the simulation to be r 0 = 1 kpc, C S0 = 18 km s −1 , and t 0 = r 0 /C S0 = 5.6×10 7 yr, respectively. The unit temperature is T 0 = µC 2 S0 /(γR) = 2×10 4 K, where µ and R are the mean molecular weight and gas constant, respectively. The normalization units are summarized in table 1.
We also introduce a non-dimensional parameter ε = γ(V 2 K /C 2 S0 ), where V K = (GM/r 0 ) 1/2 4 is the Keplerian velocity at radius r 0 . In this simulation, we adopt ε = 130 and V K = 200 km s −1 . For this parameter, the sound speed in the low-temperature region is about 5.6 km s −1 whose value is close to the observed velocity dispersion (5 − 9 km s −1 ) of the main gas components (Boulares & Cox 1990). However, this sound speed is larger than that of the molecular gas.

Initial State
The initial state is assumed to be in magnetohydrostatic equilibrium. The gas layer is initially composed of three layers: the cool equatorial layer (T = T c , |z| < z 1 ), the warm (midtemperature) layer (T = T m , z 1 ≤ |z| ≤ z 2 ), and the hot galactic halo (T = T h , |z| > z 2 ). The initial distribution of temperature is assumed to be In this study, we take T c = 0.1T 0 , T m = T 0 , T h = 10T 0 , w t = 0.015r 0 , z 1 = 0.12r 0 and z 2 = 0.24r 0 .
Although these values are not realistic for the Galactic gas disk, they would be acceptable for our first attempt to study the effects of the cold component of the galactic gas and the mechanism of the formation of dense loop-like structures. Numerical results do not depend much on the temperature of the hot galactic halo (Kamaya et al. 1997). Further discussion on the dependence of numerical results on the disk temperature will be given in section 4.5. We assume that the magnetic field is initially parallel to the equatorial plane; B = (B x (z),0,0), and is localized in the cool equatorial region (|z| < z 1 ). The magnetic field strength is determined by introducing the plasma beta (the ratio of gas pressure to magnetic pressure, hereafter denoted by β) as where 1 Here, β 0 is the plasma β at the galactic plane, and z f is the half thickness of the magnetic flux sheet. The initial density and gas pressure distributions are calculated numerically by solving the equation of magnetohydrostatic equilibrium: We have studied the evolution for five models, whose parameters are summarized in Table  2. The distributions of the initial temperature T , density ρ, gas pressure P , and magnetic pressure B 2 x /(8π) are shown in Figure 2a. Figure 2b shows the profile of the initial local pressure scale height Λ(z) = C S (z) 2 /(γg(z)) and the pressure scale height including the magnetic field 5 Λ B (z) = (1 + β(z) −1 )Λ(z) for model B (β 0 = 1), which we call the fiducial model hereafter.

Stability of the Equilibrium Model
The equilibrium state we described in section 2.2 is unstable against the Parker instability. The Parker instability in an isothermal gas has the linear maximum growth rate at a finite wavelength λ Parker ∼ 10Λ − 20Λ, where Λ is the scale height (Parker 1966). This is because small-wavelength modes are stabilized by the magnetic tension force. We analyzed the linear stability of the initial model with a normal mode method similar to that of Horiuchi et al. (1988). The linearized equations are the same as those in Horiuchi et al. (1988) and Kamaya et al. (1997). We consider the growth of a small perturbation that has a functional form δW ∝ exp(iωt + ik x x + ik y y), where W is a physical quantity (ρ, P, v z , B x ) and δW is its perturbation. The eigenvalues (ω) and the eigenfunctions are calculated numerically. Figure 3a shows the linear growth rate (iω) of the fundamental mode ) of the Parker instability as a function of the horizontal wave number k x for five cases, β 0 = 0.5,1,2,4 and 10 when z f /r 0 = 0.08 and k y = 0. Figure 3b shows that the maximum growth rate is inversely proportional to the square root of β (iω ∝ β −1/2 0 ). This result is consistent with that of the linear analysis of the Parker instability in uniform gravitational fields (Parker 1966). The most unstable wavelength (λ max ) is λ max ≃ (π/8)r 0 for model B (β 0 = 1). It is noted from Figure 2b and Figure 3 that λ max is nearly 10 times the local pressure scale height of the mid-temperature region in 0.12 ≤ z/r 0 ≤ 0.24 (λ max ∼ 10Λ ∼ 400 pc) and is much larger than that of the low-temperature cool region around z/r 0 ∼ 0.1.

Boundary Conditions and Numerical Method
We assumed free boundaries at z = Z min and z = Z max such that waves transmit freely by setting z-derivatives of all the variables vanish, and imposed periodic boundary conditions at x = X min and x = X max .
In order to trigger the Parker instability, small velocity perturbations of the form are given initially within a finite horizontal domain (|x| ≤ λ x /2) on the magnetic flux sheet, where the perturbation wavelength λ x is close to that of the most unstable wavelength for the Parker instability. Here A (= 0.01) is the maximum value of v z /C S0 in the initial perturbation. The numerical scheme we used is the Rational CIP (Cubic interpolated profile) method (Yabe & Aoki 1991;Xiao, Yabe & Ito 1996) combined with the MOC-CT method (Evans & Hawley 1988;Stone & Norman 1992). The magnetic induction equation was solved by the MOC-CT and the other equations were solved by the CIP (Kudoh, Matsumoto & Shibata 1998, 1999. The size of the simulation box is (X max − X min , Z max − Z min ) = (4r 0 , 4r 0 ) = (19H, 19H), where H is the scale height at the point where the gravity is maximum (H = C 2 S /(γg max ) = 6 0.21r 0 ). The grid sizes are ∆x = 2.5 × 10 −3 r 0 for |x| ≤ 0.8r 0 , ∆z = 2.5 × 10 −3 r 0 for |z| ≤ 0.8r 0 and they slowly increase for |x| > 0.8r 0 , or |z| > 0.8r 0 by an increment of 5% at each grid (e.g, |∆x i+1 | = 1.05|∆x i |). The number of grid points is (N x , N z ) = (472, 472). Figure 4 shows the time evolution of the fiducial model. Solid curves depict magnetic field lines. Color shows the density distribution. Arrows show velocity vectors. The overall evolution agrees with that of Matsumoto et al. (1988). That is to say, as the instability grows, the magnetic field lines bend across the equatorial plane. As the gas slides down along the undulating magnetic field lines, the rarefied regions buoyantly rise, and form magnetic loops in the later (nonlinear) phase. In the valleys of the magnetic loops, dense spur-like structures are created almost perpendicular to the galactic plane. At the top of the emerging loops, dense shell-like structures are formed (Figure 5a). These structures were not recognized in previous simulations (e.g., Matsumoto et al. 1988) partly because they assumed isothermal atmosphere without a steep density gradient at the disk-halo interface, and partly because of the lower numerical resolution. Since the downflow speed exceeds the local sound speed, strong shock waves are formed at the magnetic loop footpoints. Figure 5b shows that the shock waves heat the cool gas around z/r 0 ∼ 0.1.

Numerical Results
At the final phase (t/t 0 = 1.6), the expansion of the magnetic loops is stalled because the driving forces diminish at the top of the loop. Figure 6 shows the vertical distribution of the magnetic pressure, gas pressure, β, and density at the midpoint of the emerging loop (x/r 0 = 0.25). From this figure, it is clear that the magnetic pressure and gas pressure gradients nearly balance at the top of the emerging loop. This magneto-hydrostatic state is attained because the Parker instability is stabilized in the hot layer where the local pressure scale height becomes larger than the length of the magnetic loop. Figure 7 shows the snapshots of the density distribution for all models (β 0 = 0.5, 1, 2, 4 and 10) at the stage when the top of the magnetic loops enters the hot region (z/r 0 > 0.24). The time scale for the loop emergence is shorter for lower β, consistent with the results of the linear analysis presented in section 2.3. Figure 8 shows the vertical distributions of the density at the midpoint of the emerging loop. The density increases with height at the top of the magnetic loops in all models. In the model with β 0 = 10, the emerging loop stalls at a lower height (z/r 0 < 0.3) because the released magnetic energy is small.

Formation of Dense Loop Structure
Let us discuss why the emerging loop becomes top-heavy. Figure 9 shows the magnetic field lines at t/t 0 = 1.4 (solid curves) and t/t 0 = 0.9 (broken curves) for the fiducial model (model B). These field lines are iso-contours of the y-component of the vector potential. The magnetic field lines comoving with the plasma in ideal MHD can be identified by the value of the vector potential. Since the frozen-in condition is assumed and the numerical diffusion is negligibly small, the field lines at t/t 0 = 0.9 (broken curves) have moved with gas to the corresponding field lines at t/t 0 = 1.4 (solid curves). As the magnetic loops rise, they can accumulate the gas above the loop when the loop top is flat. Figure 9 shows that the loop shape at t/t 0 = 0.9 is favorable for the mass accumulation around the loop top. At t/t 0 = 1.4, since the curvature increases, the mass accumulated around the loop top slides down along the magnetic field lines. Flat top loops are formed when the loop length is much longer than the local pressure scale height. In our model atmosphere, since the local pressure scale height decreases with height in 0 < z/r 0 < 0.1 and increases with height in 0.1 < z/r 0 < 0.25 (see Figure 2b), the loop top tends to be flat when its height is z/r 0 ∼ 0.1. The dense shell around the loop top is formed when 0.1 < z/r 0 < 0.2, and the mass drains as the loop rises.
In order to quantitatively evaluate the density enhancement at the top of the emerging loop, we adopt the same method of analysis as Isobe et al. (2006). The velocity vector v can be divided into two components: where are the velocity components perpendicular and parallel to the magnetic field, respectively. Figure 10 shows the distribution of v and ∇ · v (Figure 10a,b), v ⊥ and ∇ · v ⊥ ( Figure  10c,d), and v and ∇ · v (Figure 10e,f) at t/t 0 = 1.2 and t/t 0 = 1.3 for the fiducial model. The distributions of ∇ · v and ∇ · v ⊥ show that the density inside the emerging magnetic loops keeps decreasing. However, ∇ · v and ∇ · v ⊥ at the top of the emerging loop are negative. Therefore, at least for the present parameters, it can be concluded that the density increases in the top region of the emerging loop. Figure 11 shows the distribution of the density near the loop top at t/t 0 = 1.1 and t/t 0 = 1.2 for the fiducial model. The clump of gas near the loop top survived due to the small curvature in the low-temperature layer and is later compressed in the high-temperature layer. Fukui et al. (2006) found two dense gas features having a loop-like shape with a length of several hundred pc within ∼ 1 kpc from the Galactic center. Figure 3a shows that the most unstable wavelength is likely to be 350 -500 pc. Our numerical simulation produced the magnetic loops whose wavelength is ∼ 400 pc and their height is ∼ 350 pc. In the observed molecular loops, the line-of-sight velocity along the loop changes linearly with the arc-length and has large gradients (∼ 30 km s −1 ). Figure 12 shows the velocity components along the outermost magnetic field line at t/t 0 = 1.4 for the fiducial model. The gradient of the downflow velocity given in Figure 12c is ∼ 20 − 30 km s −1 . Since the downflow speed exceeds the local sound speed, strong shock waves are formed at the magnetic loop footpoints. Such shocks can be the origin of the observed large velocity dispersions near the footpoints of Galactic center molecular loops. Fukui et al. (2006) estimated that the kinetic energy of the molecular loop is ∼ 10 51 erg. In our simulation, the kinetic energy flux, F k , carried by the downflow is

Estimation of the Kinetic Energy of the Emerging Loop
where ρ df is assumed to be 10 −24 g cm −3 and v df is found to be ∼ 30 km s −1 according to the numerical result for the fiducial model. The numerical results indicate that the loop width in z-direction is about 30 pc (∼ 10 20 cm) and the downflow continues for about 10 15 s ∼ 3 × 10 7 yr. When the loop thickness in y-direction is 100 pc, the total kinetic energy of the downflow toward the both footpoints is ∼ 10 51 erg. This energy is comparable to that estimated from the observation.

Effects of Cooling and Rotation
In this study, we assumed that the gas layer at the initial state is composed of the cool equatorial layer, warm layer, and hot layers. We assumed adiabatic gas but in the interstellar gas, cooling and heating play essential roles in the formation of molecular clouds (Field 1965). Kosiński & Hanasz (2006, 2007 investigated the Parker instability coupled with thermal processes (cooling and heating). They found that the Parker instability can trigger the thermal instability which form dense clouds in the valleys of the magnetic loops. In subsequent papers, we would like to include the gas cooling and heating effects.
We should consider the effect of the rotation of the disk to construct a more realistic model. In rotating disks, the Parker instability is slightly suppressed by Coriolis forces (Chou et al. 1997). Hanasz et al. (2002) presented the results of resistive 3D MHD simulations of a local part of the disk including the contribution of the Coriolis force. They demonstrated that the Parker instability, the twisting of the loops, and magnetic reconnection lead to the formation of helically twisted magnetic flux tubes. It is worth noting other 3D effects such as the growth of the interchange mode (e.g., Nozawa 2005). Isobe et al. (2006) showed that dense shells at the top of magnetic loops fragment into filaments due to the growth of the Rayleigh-Taylor instability. In order to include this effect, we need to carry out 3D MHD simulations.
In Galactic gas disks, magneto-rotational instability (MRI; Balbus & Hawley 1991) grows and drives magnetic turbulence inside the disk. In the nonlinear stage, MRI drives magnetic turbulence inside the disk. The effects of stochastic magnetic fields on the growth of the Parker instability was reported by Parker & Jokipii (2000). The results of global 3D MHD simulations of the galactic center gas disk are reported by Machida et al. (2009).

Dependence on the initial disk temperature
We studied the dependence of the numerical results on the initial disk temperature. Figure 13 shows the density distribution at the stage when the top of the magnetic loops reaches Z/r 0 ∼ 0.35 for models with lower initial disk temperature (T c = 0.05T 0 ) and higher initial disk temperature (T c = 0.2T 0 ). Other paremeters are the same as the fiducial model (model B). Numerical results indicate that the equatorial dense region at this stage is thin for the cool disk (T c = 0.05T 0 ) and thick for the warm disk (T c = 0.2T 0 ). Although the length of the loops near the equator is smaller for cooler disk, the length of the loops emerging in the halo is almost the same (∼ 400 pc). This is because the most unstable wavelength of the Parker instability is determined by the local pressure scale height.

Other Effects
Finally, we briefly mention about other effects such as cosmic rays, self-gravity and the presence of a spiral arm of the Galactic gas disk. In Galactic disks, since the cosmic ray pressure is comparable to the gas pressure, cosmic rays enhance the growth rate of the Parker instability (e.g., Parker 1966, Hanasz & Lesch 2000, 2003Kuwabara, Nakamura & Ko 2004). The effects of the cosmic rays may be more important in the Galactic center where strong activities can produce high energy particles.
In this paper, we neglected the self-gravity of the gas. When the surface density of the gas disk is large enough, Parker-Jeans instability grows (e.g., Elmegreen 1982;Nakamura, Hanawa & Nakano 1991;Kim, Ostriker & Stone 2002;Lee et al. 2004). The Jeans instability has a larger growth rate than does the Parker instability when either the magnetic field is weak or the wavelength of the perturbation is long. Franco et al. (2002) showed by 3D MHD simulations that the Parker instability creates massive clouds inside the spiral arm of the Galactic gas disk and that the dense gas accumulated around the equatorial plane forms corrugated structure. The distribution of H I gas and dust below the molecular loops found in the Galactic center (Torii et al. 2009) is consistent with the simulation by Franco et al. (2002). The density distribution in the nonlinear stage of our simulation also shows such corrugated structures.

Summary
In this paper, we have carried out 2D MHD simulations of the Galactic center gas disk consisting of the low-temperature layer and the mid-temperature layer with the overlying hot halo. We found that the numerical results reproduce basic features in the Galactic-center molecular loops as observed with NANTEN, such as the loop length, the velocity gradient along the loops, and large velocity dispersions at the footpoints of the loops.
We also discussed the reason why the top of the emerging loop becomes over-dense. This is because the effective gravity along the magnetic field lines decreases when the curvature radius of the magnetic field lines increases in the low-temperature layer. The gas survived near the loop top due to its small curvature is compressed when the loop top enters the high-temperature layer.
Finally, we suggest that the Galactic center molecular loops are analogous to the arch filaments in the solar chromosphere. The molecular loops emerge from the low temperature layer just like the dark filaments observed in Hα images of the emerging flux region of the sun. Velocity v C S0 18 km s −1       One-dimensional (z-)distribution of the logarithmic magnetic pressure P m (solid curve), gas pressure P g (broken curve), plasma beta β (dotted curve) and density ρ (dash dotted curve) at t/t 0 = 1.6 for the fiducial model (model B).