# MHD Simulations of Parker Instability Undergoing Cosmic-Ray Diffusion

## Abstract

Parker instability arises from the presence of magnetic fields in a plasma in a gravitational field such as the interstellar medium (ISM), wherein the magnetic buoyant pressure expels the gas and causes the gas to move along the field lines. The subsequent gravitational collapse of the plasma gas is thought to be responsible for the formation of giant molecular clouds in the Galaxy. The process of mixing in the ISM near the Galactic plane is investigated. The initial ISM is assumed to consist of two fluids: plasma gas and cosmic-ray particles, in hydrostatic equilibrium, coupled with a uniform, azimuthally-aligned magnetic field. The evolution of the instability is explored in two models: an isothermal exponential-declining density model and a two-layered, hyperbolic tangent temperature model. After a small perturbation, the unstable gas aggregates at the bottom of the magnetic loops and forms dense blobs. The growth rate of the instability decreases as the coupling between the cosmic rays and the plasma becomes stronger (meaning a smaller CR diffusion coefficient). The mixing is enhanced by the cosmic-ray diffusion, while the shape of the condensed gas depends sensitively on the initial equilibrium conditions. The hyperbolic tangent temperature model produces a more concentrated and round shape of clumps at the foot points of rising magnetic arches, like the observed giant molecular cloud, whereas the exponential density model gives rise to a filamentary morphology of the clumpy structure. When considering a minimum perpendicular or cross field diffusion of cosmic rays, which is often substantially smaller than the parallel coefficient , around of , the flow speed is significantly increased such that the magnetic loops extend to a greater altitude. We speculate that the galactic wind flow perpendicular to the galactic disk may be facilitated by Parker instabilities through the cross field diffusion of cosmic rays.

## 1 Introduction

Parker (1966) first identified the instability of a gravitationally stratified gaseous disk in a magneto-hydrostatic equilibrium, like the Galaxy, in response to perturbations due to magnetic buoyancy, for perturbations that occur in the magnetic field lines that lie parallel to the disk plane. This phenomenon is referred to as magnetic buoyancy instability (Hughes & Proctor 1988; Tajima & Shibata 1997) or Parker instability (Parker 1966; 1979). Solar magnetic activities such as sunspots, which is caused by the emergence of magnetic flux tubes from the interior of the sun into the solar atmosphere (Zwaan 1985, 1987), reflect this instability. The mushroom shape of the hydrogen cloud GW 123.4-1.5 (Baek, Kudoh & Tomisaka, 2008) also suggested magnetic flotation and, hence, Parker instability. In addition, giant dense CO molecular loops close to the Galactic center has been ascribed to the instability (Fukui et al., 2006). As is expected, gas aggregating at the foot points of the rising magnetic loops eventually collapses, becoming a giant cloud and a region that hosts stellar formations. On the other hand, for models galactic dynamos, Hanasz et al. (2004) found that if the effect of cosmic rays is incorporated, the efficiency of sustaining the galactic magnetic field can be improved. Otmianowska-Mazur et al. (2009) suggested that the extended, X-shaped magnetic halo structures observed in some edge-on galaxies could be attributed to the cosmic-ray driven dynamo. Actually, based on his instability, Parker (1992) introduced the idea of cosmic-ray driven dynamo to maintain the Galactic magnetic field, because the cosmic-ray pressure functions just like the magnetic pressure, capable of overpowering the gas pressure inside the magnetic flux tube and making it easier for the gas to rise.

Cosmic rays are a major component in the interstellar medium (ISM) in galaxies, whose energy density is comparable to the kinetic and magnetic energy densities of thermal plasma gas (Ferriére 2001). As high energy particles endowed with an energy of up to , of which, 90% are protons, as long as their gyroradius is significantly smaller than the characteristic spatial scales of the magnetic field, cosmic-ray particles only propagate along the magnetic field lines. When the gyroradius is comparable to the scale of magnetic field variation, cosmic rays interact strongly with the field through gyroresonant scattering by the magnetic irregularities. Although the velocity of cosmic-ray particles is close to the speed of light, the bulk motion of cosmic rays is diffusive and the bulk speed is of the order of Alfvén speed (Cesarsky 1980).

The morphology of Parker instability exhibits two modes, i.e. the undular mode and the interchange mode. The undular mode, also called Parker instability, is excited by perturbations along the magnetic field lines (the wavenumber vector of the perturbation parallel to the magnetic field ), where the falling gas creates a magnetic buoyancy greater than the restoring magnetic tension. The interchange mode, known as flute instability or magnetic Rayleigh-Taylor instability (Kruskal & Schwarzchild 1954), occurs for shorter-wavelength perturbations with perpendicular to , capable of causing two straight flux tubes to interchange and ultimately reducing the potential energy in the system. The linear growth rate of the interchange mode generally exceeds that of the undular mode owing to a rapid growth of short wavelengths. However, in the nonlinear stage, the undular mode often dominates (Matsumoto et al, 1993; Tajima & Shibata 1997). Thus, the undular mode (Parker instability) is more important than the interchange mode in astrophysical problems.

Baierlein (1983) and Matsumoto et al. (1988) performed the first one- and two-dimensional (2D), pure MHD simulations of Parker instability, respectively. In the nonlinear stage, the gas condensates to form giant clouds; in addition, a shock wave appears in the flow along the rising magnetic loop. By applying the simulation results of Matsumoto et al. (1988) to the solar atmosphere, Shibata et al. (1989b; 1990a) demonstrated that the emerging magnetic loop still expands self-similarly during the nonlinar evolution in two dimensions. Kamaya et al. (1996) adopted the supernova explosion as a perturbation in the ISM to trigger nonlinear instability. Earlier, Nozawa (1992) examined the instability deeper in the convectively unstable layer of the solar atmosphere; when considering how magnetic shear affects the flow (Hanawa et al. 1992, Nozawa 2005), although the interchange mode is stabilized, a large unstable thin structure may still arise. Shantanu et al. (1997) and Kim et al. (2000) made a further application of Parker instability to the Galactic disk without cosmic rays.

Parker’s original analysis showed that the instability has a maximum growth rate for non-zero , although non-zero is the main cause for the instability. In the 3D simulations for both solar and galactic problems given in Matsumoto & Shibata (1992) and Matsumoto et al. (1993), previous 2D results of cloud formation, presence of shock wave, and self-similar evolution are confirmed. However, the spatial and temporal scale in these studies depends on the . If a larger is applied, the magnetic loop tends to have a thinner structure and a horizontal expansion, which would suppress the upward expansion (Parker 1979). Similarly thin slice structures have been found in other 3D simulations of Kim et al. (1998, 2001, 2001) and Hanasz et al. (2002).

This study describes the evolution of Parker instability undergoing cosmic-ray diffusion using 2D simulations. Effects of the adiabatic index and the initial hydrostatic condition other than an isothermal temperature and a uniform density profile are examined. Two initial profiles are used, i.e., the hyperbolic tangent temperature model, which is often used in the solar atmosphere, and the exponential density model, which is appropriate for Galactic problems. Various perturbation results are explored. The physical parameters and initial conditions invoked are appropriate for a galactic ISM in the solar neighborhood. Section 2–6 details the governing equations, numerical algorithms, normalization, initial conditions, and grid setup. Section 7 presents the results. Section 8 draws the conclusion.

## 2 Governing Equations

We investigate the Parker instability with respect to how cosmic rays affect the fluid element in a uniformly-rotating galactic disk, under the influence of external gravity from the galactic center. A local rectangular domain representing the corotating sheet of the disk in the vertical plane is used for the simulations. The horizontal component of the radially inward external gravity is balanced by the centrifugal force, and hence, in the momentum equation, only the vertical component of external gravity and the Coriolis force are present to account for the rotation. Self-gravity from the plasma gas is not included.

A hydrodynamic approach is adopted, in which the cosmic rays and the thermal plasma are two separate fluid components of a plasma system. The plasma fluid has a mass density of , a thermal pressure of , and a cosmic-ray pressure of , all of which are threaded by a frozen-in magnetic field . The cosmic-ray energy equation is described based on the diffusion-convection equation (Drury & Völk 1981; Jones & Kang 1990, Ko 1992), which treats the cosmic rays as a hot massless fluid and neglects the momentum spectrum of cosmic rays to simplify the governing equations. The artificial separation of the cosmic rays from the plasma helps to distinguish the role played by high- and low-energy components.

The cosmic-ray diffusion-convection equation supplements the standard set of ideal MHD equations. The governing equations in the complete vector form are written as

(1) | |||||

(2) | |||||

(3) | |||||

(4) | |||||

(5) | |||||

where denotes the plasma fluid velocity; I denotes a unit tensor; denotes the adiabatic index, i.e. ratio of the heat capacity at a constant pressure to that at a constant volume, of the thermal plasma gas; refers to the adiabatic index of the cosmic rays; represents the cosmic-ray diffusivity; represents the rotation angular frequency; and is the external gravitational acceleration. Deriving this equation set involves use of the distribution function in Vlasov equation (Skilling 1975a, b, c). Such governing equations have also been used in the 2D hyperbolic tangent temperature model of Kuwabara, Nakamura & Ko (2004), and the 3D model of Lo, Ko & Wang (2011).

In addition to balancing the energy equations, the term cosmic-ray pressure affects the momentum equation Eq. (2). However, particles of cosmic rays do not interact with plasma directly; they interact with plasma via the magnetic field. On a microscopic scale, resonant scattering of Alfvén waves keeps the cosmic rays nearly isotropically distributed everywhere with respect to the thermal plasma background. A situation in which the cosmic-ray pressure possesses a gradient influences the motion of the plasma gas. Thus, in formulating the equations, the interaction of cosmic-ray particles with a thermal plasma is represented by the cosmic-ray pressure and its gradient. The transport of the cosmic-ray pressure is described by a macroscopic diffusion coefficient, , an energy-weighted mean diffusion tensor, defined as

(6) |

where denotes a unit vector along the direction.

The limit of this model is mainly that one must assume a priori knowledge of the cosmic-ray pressure and energy density that satisfies the adiabatic index . Here, Eqs. (1)–(5) are solved numerically in a local reference frame in Cartesian coordinates, whose center lies at a galactocentric radius and orbits the galaxy with a fixed angular velocity . In this local frame, the radial, azimuthal, and vertical spatial coordinates are related to the Cartesian coordinates such that , , and . In two dimensions a slice of flow is simulated, and so the rotational effect is not present.

## 3 Numerical Algorithms

A finite difference method with operator splitting is employed to solve the governing equations. Eqs. (1)–(4) are in the flux-conservative form and solved by 2-Steps Lax-Wendroff explicit method. The cosmic-ray energy equation Eq. (5) is divided into the convection and diffusion parts. The convection part is first converted to a conservative form and then also solved by 2-Steps Lax-Wendroff explicit method. The diffusion part is solved by the Bi-Conjugate Gradients Stabilized (BICGStab) implicit method.

In the 2-step Lax-Wendroff method, Eqs. (1)–(4) and the convection part of Eq. (5) are rewritten in the conservative form:

(7) |

where can be density , momentums , energy and magnetic field ; ,, are the flux of in -,-,- direction; and is the source term. The component of Eq. (7) is split into two steps:

(8) |

where superscript denotes time advection and subscript represents space grid.

The diffusion part of Eq. (5) is written as:

whose finite-difference form is expressed as , where the submatrix in and the vector is the function of time step , the cosmic-ray diffusion coefficient , the magnetic field and the cosmic-ray energy at the th iteration. The vector is the unknown:

(9) |

In order to solve in Eq. (9), the BICGStab method is employed to handle asymmetric linear systems and reduce the operations per iteration to , where is the number of unknowns in the discretized domain. This method is more efficient than the direct solution methods such as decomposition, which require operations. The BICGStab is an iteration method that uses an initial guess to find a corresponding residual , and then iterate to the -th step to an accepted value by means of bi-conjugate matrix-vector operation. For comparison, when solving the nonlinear system with Newton’s method:

the correction is determined by solving the gradient of the function, whereas in the BICGStab method,

is used for the iteration. The BICGStab treatment for cosmic-ray energy diffusion equation is an implicit method, and thus the Courant condition is not affected.

## 4 Normalization and Initial Values

The parameters and variables used in the governing equations withstand an extremely large contrast when expressed in dimensional units, possibly incurring significant errors due to the presence of ill-conditioned matrices and thus infeasible for numerical calculations. To overcome such numerical difficulties, the above equations are normalized to non-dimensional values that are close to unity. The non-dimensional values is denoted by the superscript ‘*’ while the scaling factors are denoted by a ‘0’ subscript. Some of the scaling factors also represent the quantities in equilibrium at the midplane of the Galactic disk.

The length variables are scaled based on the pressure scale height : , , and . is determined by integrating the hydrostatic equilibrium equation:

(10) |

(11) |

given the ideal gas law, the values of sound speed , the gravitational acceleration in the -direction, the ratio of plasma pressure to magnetic pressure , the ratio of plasma pressure to cosmic-ray pressure ; and . The gas at the Galactic disk plane is nearly isothermal such that . The sound speed at the disk plane is . Thus, for and , . Given , the velocity is normalized to , and the time scaling factor is defined as , i.e. The ISM density is used to normalize the density, with being about one atom per cubic cm in the ISM, or .

With the above scaling parameters, the mass equation Eq. (1) becomes

which can be converted into . Similarly, the factor at the left hand side of the momentum equation Eq. (2) is balanced by the same factor at the right hand side. The normalized plasma gas pressure is , where , and so the gradient of plasma pressure becomes .

Scaling of the cosmic-ray pressure also employs the gas pressure, . The scaling factor for the magnetic field, , is determined by , and for the Galaxy, .

The gravitational acceleration is scaled to , while the scaling factor for the rotating angular frequency is the reverse of time: . Using the same scaling parameters, the momentum equation Eq. (2), the induction equation Eq. (3) and energy equation Eq. (4) can all be converted to dimensionless values.

Finally, in Eq. (5), the cosmic-ray diffusion coefficient is normalized to , where . Theoretical calculations using path length distribution, life-time of radioactive secondary CR, etc. suggested that the diffusion coefficient for CR above 1 GeV in the Milky Way Galaxy is (Berezinskii et al. 1990; Ptuskin 2001). Due to numerical difficulties with strong magnetic field and large diffusion coefficient, this work adopts a magnetic field strength of several , like those in the solar vicinity, and a times lower than suggested. Crocker et al. (2010) suggest a lower limit of â¼ 50 on 400 pc scales near the Galactic center. However, measurements of the amplitude of the magnetic field depend on the spatial scale and the energy equipartition or pressure equilibrium among various Galactic components, which may differ by 2 orders of magnitude. Pure MHD models with strong magnetic field have been presented in Machida et al. (2009) and Takahashi et al. (2009).

## 5 Initial Equilibrium Background

We adopt two initial equilibrium backgrounds: a temperature distribution that follows a hyperbolic tangent profile and a density distribution that follows an exponential profile. Both these backgrounds bear a declining density profile with height. They differ mainly in that the hyperbolic tangent temperature has an ascending transition zone that divides the distribution into two distinct regions, while the exponential density follows a continuous decline. Given the temperature or density profile, the profiles for the density/temperature, pressure, magnetic field and other variables can be derived from the hydrostatic condition.

In the isothermal case the criteria for Parker instability to develop is . When the rising field lines grow, the flow becomes unstable. In our cases, the growth of the instability depends on additional parameters such as the width of transition region and the height of Galactic halo (disk thickness).

### 5.1 Hyperbolic Tangent Temperature

The hyperbolic tangent temperature profile is a two-temperature, layered disk (Shibata et al. 1989a), described by

(12) |

where is the sound speed, is the disk temperature and is the halo temperature. The initial dimensionless temperature is 1.0, equivalent to . Given the ideal gas law and the gravitational acceleration , the initial density profile is solved using the hydrostatic equation Eq. (10). A constant acceleration is assumed because the grid domain is small and not far way from the disk plane. The dimensionless gravitational acceleration is set as (dimensionless units), equivalent to , or , for . Acceleration can vary widely from galaxy to galaxy. Although a little higher, this value is comparable to the value derived from observations of the spatial density distribution and the velocity-distance correlation pc above the Galactic disk midplane, (Oort 1960; Bahcall 1984; Kuijken & Gilmore 1989). The initial density at the disk plane is (dimensionless units), equivalent to .

After obtaining the density profile, the plasma pressure profile is derived. The magnetic field is assumed to align in the -direction and vary with height , and . The initial dimensionless pressure at the disk is (). The initial magnetic field at the disk for is (dimensionless units), or (), This field strength is very close to the radio synchrotron measurement of 6 averaged over a radius of about 1 kpc around the Sun (Beck 2009).

The initial cosmic-ray pressure at disk for is (dimensionless units), or . The rotating angular frequency (not used in this work) is a free parameter; gives , or approximately times the angular frequency at our Sun, . The initial equilibrium background is assumed quiescent. No initial -velocity and no rotation effects are present in the system.

### 5.2 Exponential density

The exponential density case produces a nearly isothermal background that is similar to the isothermal model of Kim et al. (1998). However, since imposes singularity in the energy equation, is adopted. The exponential density profile is defined as:

(13) |

where (dimensionless units), equivalent to ; the density scale (dimensionless units) is equivalent to . The plasma pressure follows . Using Eq. (10), the temperature becomes:

(14) |

which is a constant. A minimum value of the adiabatic index is adopted in this model. Similar to the hyperbolic tangent model, the initial magnetic field and cosmic-ray pressure are defined through the parameter and , and .

## 6 Grid Setup and Perturbation

2D simulations are performed in a rectangular domain in the plane. To excite the instability, velocity perturbations are added to the quiescent background. Two perturbation forms are examined: eigenmode sinusoidal wave and random seed. The strong explosive perturbation in CR pressure developed from supernova explosions has been considered in Kuwabara et al. (2004). Notably, although supernova shocks has been believed as efficient cosmic ray accelerator, channeling at least 10% of the supernova explosion energy into cosmic rays, signature of cosmic-ray protons, the key observational evidence in support of the claim that isolated supernova remnants (SNRs) are the main accelerators of galactic CRs, has remained elusive after decades of observational effort (Allen, Houck, & Sturner 2009). Calculations of the maximum energies of shock-accelerated electrons (Reynolds & Keohane 1999) and hydrodynamic simulations of Rayleigh-Taylor instabilities in SNRs (Ferrand et al. 2010, Fraschetti et al. 2010, Wang 2011) both indicate that young SNRs could not be responsible for the highest-energy Galactic cosmic rays, unless an unrealistically high injection rate of cosmic rays, and thus an enormous energy loss in SNRs, is invoked. Precise measurements of cosmic-ray proton and helium spectra by Adriani et al. (2011) furher excludes a single power law energy distribution for the two species and thus challenged the supernova acceleration paradigm.

The sinusoidal perturbation employs a perturbing velocity described by

where is the most unstable wavelength derived from the linear analysis for (Kuwabara et al. 2004). is applied within the region of and , where (dimensionless unit). The grid consists of zones in a rectangular region of . The horizontal length of each grid zone is , while the vertical length is nonuniform:

is increased above by a ratio of 1.05, until . The upper grid domain accommodates the unperturbed magnetic field to avoid outflows across the grid boundary and spurious results.

For the random velocity perturbation (Shore & Larosa 1999), a series of random velocities and with a maximum amplitude of is added horizontally into the background (Chou et al. 2000). The grid covers an area of using zones, with and

## 7 Results

### 7.1 Pure MHD Models

We examine how cosmic-ray diffusion affects Parker instability by comparing a series of pure MHD simulations with simulations that include the cosmic-ray effect. Figure 1 illustrates those results using sinusoidal eigenmode perturbations. Figure 1(a) presents the hyperbolic tangent temperature model, while Figure 1(b) displays the exponential density model at the same epoch. In these figures a grid domain of that just lies within the maximum height of the magnetic loops is shown. Logarithmic values of variables in a dimensionless unit are presented in color contours. White solid lines represent the magnetic field lines, and small white arrows refer to the velocity vectors. The superscript ’’ denoting dimensionless values is omitted hereafter. The exponential model shows a faster growth of the instability because the decline in the gas density facilitates the instability.

### 7.2 Cosmic-Ray Diffusion Coefficient

Figure 2 compares the maximum height of the magnetic loops versus the cosmic-ray diffusion coefficient in both models at , when the numerical time step becomes very restrained. A smaller CR diffusion coefficient means stronger coupling between the cosmic rays and the plasma. A maximum of ( is examined. In both models, the instability becomes more prominent with an increasing , and the mixing extends to a similarly maximum height of . Because the gradient of cosmic-ray pressure declines with an increasing diffusivity, cosmic-ray diffusion facilitates flow instability. Kuwabara et al. (2004) studied the instability in the hyperbolic tangent temperature case. For small-, the instability is impeded by the CR pressure gradient force that interferes with the falling motion of the matter, while for large-, the magnetic loop can rise to larger scales. Our simulations reach a similar outcome. In Kuwabara et al (2004), the case of strong explosive perturbation in CR pressure presumably from supernova yielded a larger growth rate for a smaller diffusion coefficient, but such a reversed correlation only occurs in the early stage; in the later stage the growth rate becomes smaller when compared to that of a large diffusion coefficient model. In the 3D models of Lo, Ko & Wang (2011) employing sinusoidal perturbations, however, the instability growth first rises with but then drops when , because under strong diffusion, the falling of the gas becomes supersonic; as a result of shock, the cosmic rays are redistributed with a smoother pressure gradient, which tends to stabilize the flow.

### 7.3 Parameter and

The scale height is set for . Varying modifies the length scale, resulting in different flow structures. Figure 3 presents snapshots of the hyperbolic tangent temperature model with a tenfold gas pressure, ; the flow extends to , which displays more prominent mixing than the case (). Other parameters used are ; ; ; (dimensionless units, or ; ( or ; (); and .

### 7.4 Parameter and

The temperature distribution is assumed to be uniform in the exponential density model. In this case, modifying and changes only the isothermal temperature value. Notably, the parameter should also influence the scale length (see the scaling for ). However, because cosmic rays largely diffuse along the magnetic field lines, varying does not significantly affect the scale height.

The isothermal case is approximated by using a minimum adiabatic index of . In an ISM, should be significantly smaller than the ideal value and close to , because thermal instability smooths out the temperature gradient. However, because the diffusion of cosmic rays along the magnetic field lines is only slight and the instability growth expedites for more condensed gas, thermal instability is expected to be deterred when cosmic-ray diffusion is present (Shadmehri 2009). This phenomenon may resemble the effect of increasing above the typical value that is used in most of the simulations for Parker instability. According to our results, instability is dampened as increases; with , the flow remains quiescent at (Figure 4a).

In the hyperbolic tangent temperature model, an increasing width of transition region makes the instability less feasible (Figure 4b). Suppression of the instability is attributed to the rapid rise of temperature and the flattening of density near the transition region.

### 7.5 Random Perturbation

The case using random velocity perturbation may reflect the circumstances closer to the actual situation in an ISM. In this case, the development of instability is slower than the case using eigenmode perturbation, Figure 5 reveals that the exponential density model exhibits slightly higher floatation of the magnetic loops, but similarly strong mixing and filamentary structures are induced in both models.

### 7.6 Trebly Sinusoidal Perturbation with

The perpendicular or cross field lines diffusion coefficient (see Eq. 6) is often substantially smaller than the parallel coefficient , which is only around of (Ryu et al. 2003). We found that when incorporating cosmic-ray cross field line diffusion , the magnetic field lines are aligned more vertically to the disk. As such a process may be related to the acceleration of galactic wind, a new initial condition using the same equilibrium background is invoked to examine the case involving .

A grid of zones in a domain of is adopted, with each zone having a size identical to the case of random perturbation. The applied velocity perturbation is similar to our previous sinusoidal eigen mode perturbation, but within a finite rectangular region of , where , along with other similar parameters: ; ; ; ; ; ; and .

Figure 6 shows the exponential case at the dimensionless time . Three equally-sized loops arise and extend to . The flow pattern resembles the hyperbolic tangent case in Kuwabara et al. (2004). The interaction between loops becomes more pronounced at later times; the central one becomes compressed by two outer adjacent loops, while the loops continually rise to a height of . At , the flow reaches a maximum vertical height and becomes saturated. Different from the case without , the flow velocity is significantly increased such that the loops rise to a much higher altitude, . Therefore, even at a value of of , the effect of on the mixing appears significant (Figure 6). Although such a result contradicts the estimate made by linear analysis (Ryu et al. 2003), it indicates that Parker instability is likely to contribute to the galactic wind acceleration perpendicular to the galactic disk (Wiegelmann, Schindler, & Neukirch 1997).

## 8 Conclusion

This study has elucidated the development of Parker instability undergoing cosmic-ray diffusion using 2D MHD simulations. Exactly how the initial conditions affect the growth of instability is examined. Two equilibrium backgrounds are constructed based on hyperbolic tangent temperature and exponential density, with variations in the physical parameters and perturbations.

Our simulations gives rise to a structure that is closely resembles the long, thin wave-like sheared helmet plumes in the solar corona, caused by magnetic buoyancy (Shore & Larosa 1999). The perpendicular component of cosmic-ray diffusion coefficient is included as well. When is only of the parallel component , despite the small development, vertical magnetic structures and outward flow arise, particularly in the exponential density model. The substantial increases in the velocity and altitude of the unstable flow reveal the importance of cosmic rays to the galactic wind. In general, in the exponential density model, the morphology of the clumpy structure is filamentary. Conversely, in the hyperbolic tangent temperature model, a more concentrated and round shape of clumps like the giant molecular cloud are observed at the foot points of rising magnetic arches. The growth of instability in the hyperbolic tangent model is less rapid than that of the exponential density case since the pressure in the exponential density decreases faster with an increasing height.

The small gradient of cosmic ray pressure along the z direction caused by diffusion explains why an increase in facilitates the growth of the unstable undular mode. As a result, the gas falls down more rapidly, and adjacent loops join together to form a large loop-like bubble.

Exponential density with the adiabatic index (i.e. isothermal) produces the highest magnetic loops and most flow instability. A decreasing leads to more condensed gas and, ultimately, a more unstable flow. While examining a situation without cosmic-ray diffusion, Parker indicated that the criterion for instability is . Although the undular mode is expected to be suppressed when , numerical results in our cases using indicate that the high adiabatic index still produces a non-uniform density distribution before the simulations end due to stringent numerical conditions. Therefore, the instability may still occur for large if a stronger perturbation such as a supernova explosion is invoked.

The hyperbolic tangent temperature model yields a density distribution that increases with height, whose steepness depends on the width of the transition zone . For a smaller , instability is increased as the density gradient correspondingly diminishes near the unstable regions.

When a larger ratio of plasma pressure to magnetic pressure is applied, the pressure inside the flux tube diminishes, thus the instability growth is facilitated. Also, as the length scale decreases accordingly, unstable structures are observed on a smaller scale. Contrarily, a small yields a larger scale height and so complicates the mixing of an unstable flow.

Finally, the ratio of plasma pressure to cosmic-ray pressure does not influence the length scale since the diffusion of the cosmic-ray pressure does not contribute to the scale height.

## Acknowledgment

We thank the National Center for High-Performance Computing Center for providing computing facilities and the National Science Council of Taiwan for funding this project under Grants NSC 98-2112-M-033-002, NSC 98-2811-M-033-014 and NSC 98-2811-M-033-020. CMK is supported by the Taiwan National Science Council Grants NSC 98-2923-M-008-01-MY3 and NSC 99-2112-M-008-015-MY3. We are grateful to Jongsoo Kim for helpful comments.

### References

- Adriani O. et al., 2011, Science, 332, 69
- Allen G. E., Houck J. C., Sturner S. J., 2008, ApJ, 683, 773
- Baek C. H., Kudoh T., Tomisaka K., 2008, ApJ, 682, 434
- Baierlein R., 1983, MNRAS, 205, 669
- Bahcall J. N., 1984, ApJ, 276, 169
- Beck R., 2009, Astrophys. Space Sci. Trans., 5, 43
- Berezinskii V. S., Bulanov S. V., Dogiel V. A., Ginzburg V. L., Ptuskin V. S. 1990, Astrophysics of Cosmic Rays, ed. Berezinskii V. S.Ginzburg V. L.(New York: North-Holland)
- Cesarsky, C. J., 1980, ARAA, 18, 289
- Chou W., Matsumoto R., Tajima T., Umekawa M. Shibata K., 2000, ApJ, 538, 710
- Crocker R. M., Jones D. I., Melia F., Ott J., Protheroe R. J., 2010, Nature, 463, 65
- Drury L. O. C., Völk H. J., 1981, ApJ, 248, 344 Science Reviews, 99, 329
- Ferrand G., Decourchelle1 A., Ballet J., Teyssier R., Fraschetti, F., 2010, A&A, 509, L10
- Fraschetti F., Teyssier T., Ballet J., Decourchelle A., 2010, A&A, 515, 104
- Ferriére K. M., 2001, Rev. Mod. Phys., 3, 1031
- Fukui Y., et al., 2006, Science, 314, 106
- Hanasz M., Otmianowska-mazur K., Lesch H., 2002, A&A, 386, 347
- Hanasz M., Kowal G., Otmianowska-Mazur K. and Lesch H., 2004, ApJ, 605, L33
- Hanawa T., Matsumoto R., Shibata K., 1992, ApJ, 393, L71
- Hughes D. W., Proctor M. R. E., 1988, Annual Review of Fluid Mechanics, 20, 187
- Jones T. W., Kang H., 1990, ApJ, 363, 499
- Kamaya H., Mineshige S., Shibata K., Matsumoto R., 1996, ApJ, 458, L25
- Kim J., Franco J., Hong S. S., Santillan A., Martos M. A., 2000, ApJ, 531, 873
- Kim J., Hong S. S., Ryu D. Jones T. W., 1998, ApJ, 506, L139
- Kim J., Ryu D., Jones T. W., 2001, ApJ, 557, 464
- Kim W., Ostriker E. C., Stone J. M., 2002, ApJ, 581, 1080
- Kruskal K. D., Schwarzchild M., 1954, Proc. R. Soc. London, Ser A., 223, 348
- Kuijken K., Gilmore G., 1989, MNRAS, 239, 651
- Ko C. M., 1992, A&A, 259, 377
- Kuwabara T., Nakamura K., Ko, C. M., 2004, ApJ, 607, 828
- Lo Y. Y., Ko C. M., Wang C. Y., 2010, Computer Physics Communications, 182, 177
- Machida M., Matsumoto R., Nozawak S., Takahashi K., Fukui Y., Kudo N., Torii K., Yamamoto H., Fujishita M., Tomisaki K., 2009, PASJ, 61, 411
- Matsumoto R., Horiuchi T., Shibata K., Hanawa T., 1988, PASJ, 40, 171
- Matsumoto R., Shibata K., 1992, PASJ, 44, 167
- Matsumoto R., Tajima T., Chou W., K., A. Okubo, Shibata K., 1998, ApJ, 493, 414
- Matsumoto R., Tajima T., Shibata K., Kaisig M., 1993, ApJ, 414, 357
- Nozawa S., Shibata K., Matsumoto R., Sterling A. C., Tajima T., Uchida Y., Ferrari A., Rosner R., 1992, ApJS, 78, 265
- Nozawa S., 2005, PASJ, 57, 995
- Oort J. H., 1960, Bull. Astr. Inst. Netherlands, 15, 45
- Otmianowska-Mazur K., Soida M., Kulesza-Źydzik B., Hanasz M. Kowal G., 2009, ApJ, 693, 1
- Parker E. N., 1966, ApJ, 145, 881
- Parker E. N., 1979, Cosmical Magnetic Fields, Oxford: Clarendon Press
- Parker E. N., 1992, ApJ, 401, 137
- Ptuskin V. S. 2001, Space Sci. Rev., 99, 281
- Reynolds S. P., Keohane J. W., 1999, ApJ, 525, 368
- Ryu D., Kim J., Hong S. S., Jones T. W., 2003, ApJ, 589, 338
- Shantanu B., Mouschovias T. Ch., Paleologou E. V., 1997, ApJ, 480, 55
- Shibata K., Nozawa S., Matsumoto R., Sterling A. C., Tajima T., 1990, ApJ, 351, L25 (1990a)
- Shibata K., Matsumoto R., 1991, Nature, 353, 633
- Shibata K., Tajima T., Matsumoto R., Horiuchi T., Hanawa T., Rosner R., and Uchida T., 1989, ApJ, 338, 471 (1989a)
- Shibata K., Tajima T., Steinolfson R. S., Matsumoto R., 1989, ApJ, 345, 584 (1989b)
- Shore S. N., Larosa T. N., 1999, ApJ, 521, 587
- Shadmehri M., 2009, MNRAS, 397, 1521
- Skilling J., 1975, MNRAS, 172, 577
- Skilling J., 1975, MNRAS, 173, 245
- Skilling J., 1975, MNRAS, 173, 255
- Tajima T., Shibata K., 1997, Plasma Astrophysics, Reading, Massachusetts: Addison-Wesley
- Takahashi K., Nozawa S., Matsumoto R., Machida M., Fukui Y., Kudo N., Torii K., Yamamoto H., and Fujishita M., 2009, PASJ, 61, 957
- Wang, C.-Y., 2011, MNRAS, 415, 83
- Wiegelmann T., Schindler K., Neukirch T., 1997, Solar Phys., 180, 439
- Zwaan C., 1985, Sol. Phys., 100, 397
- Zwaan C., 1987, ARAA, 25, 83