# The Parker Instability in a Thick Galactic Gaseous Disk: I. Linear Stability Analysis and Nonlinear Final Equilibria

## 1 INTRODUCTION

In a series of seminal papers, Parker (1966, 1967, 1969) discussed the stability of a magnetized interstellar system with cosmic rays, and immersed in an external gravitational field. To build an equilibrium state in a plane-parallel density distribution, he assumed that i) the initial magnetic field is parallel to the galactic plane, ii) the gravitational acceleration is constant, iii) and the vertical pressure distributions for the gas, cosmic rays and magnetic field are simply described by an exponential function with the same scale height. Using a normal mode analysis he found that such a system is unstable if the adiabatic index of the gas is below a certain critical value. When the perturbation wavevectors are confined to the two-dimensional (2D) plane defined by the directions of the initial magnetic and gravitational fields, the critical value is defined by $`\gamma _{\mathrm{cr},u}=(1+\alpha +\beta )^2/(1+0.5\alpha +\beta )`$, where $`\alpha `$ is the ratio of the magnetic-to-gas pressures, and $`\beta `$ is the ratio of cosmic-ray-to-gas pressures (Parker 1966). When the wavevectors are allowed to have all three-dimensional (3D) components, the critical adiabatic index becomes $`\gamma _{\mathrm{cr},m}=1+\alpha +\beta `$ (Parker 1967).

Later, these (2D and 3D) types of perturbations were classified as the undular (2D) and interchange (Hughes & Cattaneo 1987) or mixed (3D) modes, (Matsumoto et al. 1993). Using the “energy principle” method, Lachièze-Rey et al. (1980) also found a generalized form for the critical adiabatic index (eq.  in their paper), which is basically the same one of the mixed mode. The critical adiabatic index for this mixed mode is in general smaller, and more restrictive, than that of the undular mode. Given that cooling times in the diffuse interstellar medium (ISM) are shorter than the timescales for the instability, Parker used the isothermal value for the adiabatic index, $`\gamma =1`$, and concluded that the equilibrium state of the general ISM is then unstable.

The undular instability, which promotes the formation of high-density structures, is eventually stabilized by the tension of the distorted field lines. Using mass invariance along flux tubes, Mouschovias (1974) obtained the 2D final equilibrium state of the original one-dimensional Parker model. Even when this new 2D equilibrium stage is in turn unstable against 3D perturbations (see Asséo et al. 1978; Asséo et al. 1980), the undulation pattern along the initial field lines persists in the nonlinear 3D evolution of the instability (Kim et al. 1998). Therefore, the final 2D equilibrium state can be helpful in visualizing the resulting large-scale structure of the ISM.

The three original assumptions made by Parker described above are obvious idealizations, and some of them have been modified in subsequent studies. The first assumption, a well ordered field that is parallel to the galactic plane, is not really sustained by observations (except near the midplane). The interstellar $`𝐁`$-field has a bisymmetric spiral field configuration (see Heiles 1996; Indrani & Deshpande 1998; Vallee 1998), and random components with cell sizes of the order of 50 pc (e.g. Rand & Kulkarni 1989). Also, the transition between the gaseous disk and the halo is very broad and has a complex structure with vertical field components (see Boulares & Cox 1990 for a discussion of the support provided by the tension of curved field lines). Thus, the plane-parallel field assumption is only valid as an average field configuration, but it is very difficult to relax in both analytical and numerical treatments of the problem. As a variation to the simplest plane-parallel field scheme, Hanawa, Matsumoto & Shibata (1992) derived the unstable modes in a skewed magnetic field whose direction is still horizontal, but the field direction changes with distance from the midplane (i.e., the $`x`$ and $`y`$ components of the field vary with height, but $`B_z`$ is always equal to zero). For such a field configuration, the instability tends to form structures with the scales of giant molecular clouds.

The second assumption, a constant gravity, has been relaxed in more recent studies. The galactic gravitational field varies in a nearly linear fashion near the midplane (the linear approximation is excellent for $`[z]150`$ pc; Oort 1965; Bahcall 1984; Bienaymé, Robin & Crézé 1987; Kuijken & Gilmore 1989), and two different functions, linear and tanh ($`z/H`$), have been considered by different authors (Giz & Shu 1993; Kim, Hong & Ryu 1997; Kim & Hong 1998). Since gravity is the driving force of the instability, a variation of the functional form for the acceleration has a direct impact on the properties of the unstable modes (i.e., growth rates, length scales, and parity). For the constant gravity case, the weight of a gas parcel is the same regardless of its $`z`$-position, but the acceleration is discontinuous at $`z=0`$. Thus, the flow cannot move across the midplane and the only allowed modes are those with even parity. The resulting structures are then distributed symmetrically with respect to $`z=0`$, and are called midplane symmetric (MS). In the case of the other two functions (linear and tanh), the acceleration is continuous at the midplane and the weight of the gas increases with $`[z]`$. Thus, the odd parity solutions, or midplane antisymmetric (MA) modes, also appear, and the growth times are shorter than those of the uniform gravity case (e.g. Giz & Shu 1993; Kim & Hong 1998).

The third assumption, a disk with a single gas component, has not been modified in any of the recent studies but it also has to be revised. The actual ISM structure is very complex, and has several gas components ranging from the cold molecular phase to a hot and highly ionized plasma (e.g. Kalberla & Kerp 1998). Each component, in turn, has several sub-components with their own set of representative values for midplane densities and scale heights (the velocity dispersion of most components, however, seems to be equal to $`9`$ km s<sup>-1</sup>: see Boulares & Cox 1990). For instance, the atomic hydrogen phase could be divided into one cold H I component and two warm H I components (e.g. Bloemen 1987; Boulares & Cox 1990; McKee 1990; Spitzer 1990). To further complicate the situation, the vertical distributions for the magnetic field and cosmic rays do not seem to follow the stratification of the main gas components. Many of the system properties remain largely unknown and, depending on the assumed temperature and $`𝐁`$-field distributions, the resulting magnetohydrostatic (MHS) equilibrium configurations can be either stable or unstable to the Parker instability (e.g. Bloemen 1987; Boulares & Cox 1990; Martos & Cox 1994; Franco, Santillán & Martos 1995; Kalberla & Kerp 1998). Thus, a quantitative stability analysis for this type of multi-component gaseous disk is required.

In this paper we address this issue and investigate the 2D stability of an extended, multi-component, magnetized disk with a “realistic” gravitational acceleration. We use the vertical equilibrium model for the warm magnetized system that has been discussed by Martos (1993), Martos & Cox (1994, 1998) and Santillán et al. (1999a). This equilibrium configuration is based on the observed distributions of: i) the vertical acceleration of the gravitational field in the solar neighborhood (Bienaymé, Robin & Crézé 1987), and ii) the density distributions of the gaseous components (Boulares & Cox 1990). Using the normal mode analysis, here we find that an isothermal extended disk is unstable with respect to the undular mode, and derive the resulting linear growth rates. We also derive the final 2D equilibrium state for both the MS and MA modes. The non-linear evolution is followed with the aid of 2D magnetohydrodynamic (MHD) numerical experiments and the results will be presented in the accompanying paper by Santillán et al. (1999b, hereafter Paper II).

The plan of the present paper is as follows. In §2, we describe the initial MHS equilibrium state, and perform the normal mode analysis. The dispersion relations for the unstable undular modes are then discussed. In §3, the nonlinear final equilibria of the undular modes are presented, and a summary and discussion of the results are given in §4.

## 2 Normal Mode Analysis

The linear stability analysis can be performed with either the “energy principle” method (Bernstein et al. 1958) or with the usual “normal mode” analysis. The energy principle method provides the critical adiabatic index in a relatively easy way (e.g. Zweibel & Kulsrud 1975; Asséo et al. 1978; Asséo et al. 1980; Lachièze-Rey et al. 1980), but does not allow to derive the resulting dispersion relations. Here we want to find the dispersion relations and, hence, derive them with the usual normal mode analysis.

### 2.1 Isothermal Magnetohydrodynamic Equations

The temperature distribution of the gaseous disk is largely unknown but, given that the velocity dispersion of the main gas components are similar, we identify the velocity dispersion with the gas sound speed and define our model as an isothermal disk stratification. Thus, we do not differentiate between the thermal and kinetic pressures, and both are gathered together in a single pressure term with either a constant velocity dispersion or a constant effective temperature. The particular values for the resulting model velocity dispersion and effective temperature are given below.

The dynamics of a magnetized isothermal plasma immersed in a gravitational field is described by the MHD equations,

$$\frac{\rho }{t}+(\rho 𝐯)=0,$$
(1)   
$$\rho \left(\frac{𝐯}{t}+𝐯𝐯\right)=\left(\rho a^2+\frac{B^2}{8\pi }\right)+\frac{1}{4\pi }𝐁𝐁+\rho 𝐠,$$
(2)   
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁),$$
(3)   
where $`a`$ (= constant) is the isothermal sound speed, and the rest of the symbols have their usual meanings. We use a Cartesian coordinate system ($`x,y,z`$), whose axes are defined parallel to the radial, azimuthal, and vertical directions, respectively. We perform the analysis in the $`yz`$ plane and assume that the gravitational acceleration has only a vertical component, $`𝐠=[0,0,g(z)]`$.

### 2.2 Initial Equilibrium Configuration

After Parker (1966) introduced the simplified exponential equilibrium model for the gaseous disk, several authors have built more complex and realistic approximations for the actual ISM structure (e.g. Badhwar & Stephens 1977; Bloemen 1987; Boulares & Cox 1990; Kalberla & Kerp 1998). These models are based on the observed stratifications for the gas, cosmic rays, and magnetic and gravitational fields in the solar neighborhood. Here we use, as the initial equilibrium state in our stability analysis, the MHS equilibrium configuration discussed originally in Martos (1993). This initial model uses the vertical distributions for the density and gravitational acceleration described by Boulares & Cox (1990) and Bienaymé, Robin & Crézé (1987), respectively.

The density stratification is

$`n_0(z)`$ $`=`$ $`0.6\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(70\text{pc})^2}}\right]+0.3\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(135\text{pc})^2}}\right]+0.07\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(135\text{pc})^2}}\right]`$ (4)   
$`+`$ $`0.1\mathrm{exp}\left[{\displaystyle \frac{|z|}{400\text{pc}}}\right]+0.03\mathrm{exp}\left[{\displaystyle \frac{|z|}{900\text{pc}}}\right]\mathrm{cm}^3.`$   
The midplane value is $`n_0(0)1.1\mathrm{cm}^3`$ and, for a plasma with 10% He, the corresponding mass gas density stratification is $`\rho _0(z)=1.27m_\mathrm{H}n_0(z)`$, where $`m_\mathrm{H}`$ is the mass of a hydrogen atom. Figure 1 shows the distribution of each gas component (representing the contributions of H<sub>2</sub>, cold H I, warm H I in clouds, warm intercloud H I, and warm diffuse H II). The molecular and cold atomic phases are the dominant ISM mass components near the midplane, whereas the warm intercloud H I and warm diffuse H II are the most important gas layers beyond $`z300`$ pc. The extended H II component was originally detected in absorption against the Galactic synchrotron background (Hoyle & Ellis 1963), and later it was reported in hydrogen recombination emission (Reynolds 1989). This ionized gas is a major component of the ISM, which has been usually ignored in previous modeling, and its surface density is about a third of that of the H I component at the solar neighborhood. The power requirements to ionize this layer are comparable to that available from supernovae.

Because of the inclusion of the extended components, with scale heights larger than 300 pc, our model is referred to as the thick gaseous disk model. The resulting effective scale height is defined by the total gas column density as

$$H_{\mathrm{eff}}=\frac{1}{n_0(0)}_0^{\mathrm{}}n_0(z)𝑑z166\mathrm{pc}.$$
(5)   
In the case of the gravitational field, the observationally derived acceleration at the solar neighborhood can be fitted by (Martos 1993)

$$g(z)=8\times 10^9\left[10.52\mathrm{exp}\left(\frac{|z|}{325\text{pc}}\right)0.48\mathrm{exp}\left(\frac{|z|}{900\text{pc}}\right)\right]\mathrm{cm}\mathrm{s}^2.$$
(6)   
This gravitational acceleration is similar to the one derived by Kuijken & Gilmore (1989), and requires less local dark matter content than the ones derived by Oort (1965) and Bahcall (1984).

Given these two basic building blocks, $`\rho _0(z)`$ and $`g(z)`$, the initial equilibrium configuration is constructed by assuming that the gas is isothermal and the initial magnetic field is parallel to the galactic plane. The effects of cosmic rays are not explicitly included here because the results may depend on the assumptions made. For instance, in contrast to the effects of the isotropic cosmic ray pressure considered by Parker (1966), Nelson (1985) showed that an anisotropic cosmic ray pressure may tend to stabilize the gas layer. For simplicity, then, we gather the non-thermal pressures into a single term represented by the magnetic pressure (i.e., we assume that the sum of the cosmic ray and magnetic pressures is contained in the magnetic term). Then, the MHS equilibrium for the gas-field-gravity system is given by

$$\frac{d}{dz}P_0(z)=\frac{d}{dz}\left[\rho _0(z)a^2+\frac{B_0^2(z)}{8\pi }\right]=\rho _0(z)g(z),$$
(7)   
where $`P_0(z)`$ is the total pressure of the system (thermal plus magnetic). This equation defines the stratification of the magnetic field. For completeness, because our system is finite, we set the additional boundary condition $`P_0(z=10\mathrm{k}\mathrm{p}\mathrm{c})=0`$, and the system pressure is computed with the integral

$$P_0(z)=_z^{10\mathrm{k}\mathrm{p}\mathrm{c}}\rho _0(z)g(z)𝑑z.$$
(8)   
Given the total pressure and the strength of the magnetic field (including both ordered and random components) at midplane, $`P_0(0)3\times 10^{12}`$ dyne cm<sup>-2</sup> and $`B_0(0)5\mu `$G (Boulares & Cox 1990; Heiles 1996), the resulting isothermal sound speed (from $`P_0(0)=1.27m_\mathrm{H}n_0(0)a^2+B_0^2(0)/8\pi `$, with $`n_0(0)=1.1`$ cm<sup>-3</sup>) is $`a=8.4`$ km s<sup>-1</sup>. Thus, the sound speed value is very similar to the observed velocity dispersion of the main gas components (within 5 to 9 km s<sup>-1</sup>; Boulares & Cox 1990), and the corresponding effective disk temperature is $`T_{\mathrm{eff}}=10900`$ K. This is called the “warm” magnetic disk model and its properties are discussed by Martos (1993), Martos & Cox (1994, 1998) and Santillán et al. (1999a). Figure 2 shows the distributions of the thermal, magnetic, and total pressures as functions of distance from the galactic plane. The maximum of the magnetic pressure is not centered at $`z=0`$ because the field stratification is derived from MHS equilibrium. This warm magnetic disk model is Parker unstable because the gas is almost entirely supported by the magnetic field above $`z200`$ pc.

There is high-latitude H I gas with velocity dispersions of 35 km s<sup>-1</sup> (Kulkarni & Fich 1985), and halo gas with up to 60 km s<sup>-1</sup> (Kalberla et al. 1998). The inclusion of these additional gas components with different velocity dispersions in our analysis is beyond the scope of the present paper but, as sketched in §4, we will address this issue in a future study.

### 2.3 Linearized Perturbation Equations

We limit the present discussion to perturbations in the $`yz`$ plane (i.e., in the plane defined by the directions of the initial magnetic and gravitational fields). Due to this limitation, only the undular modes are allowed and we follow the procedure described by Kim, Hong & Ryu (1997) and Kim & Hong (1998) to derive the properties of the instability. Given that the velocities in the initial model are equal to zero, we denote by $`𝐯`$, $`\delta \rho `$, $`\delta 𝐁`$ the infinitesimal perturbations in velocity, density, and magnetic field, respectively. The perturbed state is then described by

$$𝐯;\rho =\rho _0+\delta \rho ;𝐁=B_0\widehat{e}_y+\delta 𝐁.$$
(9)   
Inserting these perturbed variables in equations (1-3), and keeping only the first-order terms for the perturbations, the linearized perturbation equations become

$$\frac{}{t}\delta \rho +v_z\frac{d\rho _0}{dz}+\rho _0\left(\frac{v_y}{y}+\frac{v_z}{z}\right)=0,$$
(10)   
$$\rho _0\frac{v_y}{t}+\frac{}{y}(a^2\delta \rho )\frac{1}{4\pi }\frac{dB_0}{dz}\delta B_z=0,$$
(11)   
$$\rho _0\frac{v_z}{t}+\frac{}{z}(a^2\delta \rho )+\frac{1}{4\pi }\frac{dB_0}{dz}\delta B_y+\frac{1}{4\pi }B_0\frac{}{z}\delta B_y\frac{1}{4\pi }B_0\frac{}{y}\delta B_z+g\delta \rho =0,$$
(12)   
$$\frac{}{t}\delta B_y+B_0\frac{v_z}{z}+\frac{dB_0}{dz}v_z=0,$$
(13)   
$$\frac{}{t}\delta B_zB_0\frac{v_z}{y}=0.$$
(14)   
The coefficients of equations (10-14) do not depend explicitly on $`y`$ and $`t`$, and the perturbations can be Fourier-decomposed with respect to these variables

$$\left[\begin{array}{c}\delta \rho (y,z;t)\\ v_y(y,z;t)\\ v_z(y,z;t)\\ \delta B_y(y,z;t)\\ \delta B_z(y,z;t)\end{array}\right]=\left[\begin{array}{c}\delta \rho (z)\\ v_y(z)\\ v_z(z)\\ \delta B_y(z)\\ \delta B_z(z)\end{array}\right]\mathrm{exp}(i\omega tik_yy),$$
(15)   
where $`i\omega `$ is the growth rate and $`k_y`$ is the wavenumber along the $`y`$-direction. Inserting these decomposed forms into the perturbation equations (10-14), and combining them we obtain the reduced equation,

$$f\frac{d^2v_z}{dz^2}+\frac{df}{dz}\frac{dv_z}{dz}+hv_z=0,$$
(16)   
where the functions $`f`$ and $`h`$ are defined by

$$f=2(\omega ^2k_y^2a^2)\frac{B_0^2}{8\pi }+\omega ^2\rho _0a^2,$$
(17)   
$$h=(\omega ^2k_y^2a^2)\left(\omega ^2\rho _02k_y^2\frac{B_0^2}{8\pi }\right)\omega ^2\rho _0\frac{dg}{dz}+k_y^2g\frac{d}{dz}\left(\frac{B_0^2}{8\pi }\right).$$
(18)   
The factor $`dg/dz`$ appearing in the second term of $`h`$ is introduced by taking the derivative with respect to $`z`$ on both sides of the MHS equation (7), and then making the appropriate substitutions. This results in a more compact form for the $`h`$ function (and we do not need to calculate numerically a second-order derivative term). With the transformation $`\mathrm{\Psi }=v_zf^{1/2}`$, equation (16) can be rearranged to

$$\mathrm{\Psi }^{\prime \prime }+\left[\frac{1}{4}\left(\frac{f^{}}{f}\right)^2\frac{1}{2}\left(\frac{f^{\prime \prime }}{f}\right)+\frac{h}{f}\right]\mathrm{\Psi }=0,$$
(19)   
where the prime superscript () denotes the derivative with respect to $`z`$. Given the complicated functional forms for $`\rho _0(z)`$, $`g(z)`$, and $`B_0(z)`$, one cannot perform further simplications of equation (19).

The required boundary conditions (BCs) are: $`\mathrm{\Psi }=0`$ at an upper boundary $`z=z_{\mathrm{node}}`$, and $`\mathrm{\Psi }=0`$ or $`d\mathrm{\Psi }/dz=0`$ at the midplane, $`z=0`$. The first condition at the midplane, $`\mathrm{\Psi }=0`$ at $`z=0`$, generates the even parity MS solutions, whereas the second one corresponds to the odd parity MA solutions (e.g. Horiuchi et al. 1988; Giz & Shu 1993).

### 2.4 Dispersion Relations

The dispersion relations are found with the method described in the Appendix of Kim et al. (1997). The method is a numerical procedure to find, for a given wavenumber, an eigenvalue ($`i\omega `$) which satisfies the imposed BCs. Our equilibrium configuration, as stated in §2.2, turns out to be Parker unstable, and we find eigenvalues $`i\omega `$ that are real and positive. The resulting dispersion relations are shown in Figure 3 for five cases whose upper boundaries are placed at $`z`$-locations ranging from $`z=9H_{\mathrm{eff}}`$ to $`z=30H_{\mathrm{eff}}`$. The growth rates, wavenumbers, and nodal points in Figure 3 are normalized as:

$$\mathrm{\Omega }=i\omega \frac{H_{\mathrm{eff}}}{a},\nu _y=k_yH_{\mathrm{eff}},\zeta _{\mathrm{node}}=\frac{z_{\mathrm{node}}}{H_{\mathrm{eff}}}.$$
(20)   
Since gravity has small values near the midplane, the dispersion relations are not sensitive to the midplane boundary conditions. Thus, the solutions are degenerate with respect to parity, and the growth rates are nearly the same for both the MS and MA modes. The plotted dispersion relations are for the principal $`z`$-modes, whose MS nodal points are located at midplane and $`\zeta _{\mathrm{node}}`$. The nodal points of the MA modes are located at $`\zeta _{\mathrm{node}}`$ and $`\zeta _{\mathrm{node}}`$. For the lower nodal point in Figure 3, $`\zeta _{\mathrm{node}}=9`$ ($`z_{\mathrm{node}}=1.5`$ kpc), the fastest growth time is about $`6.2\times 10^7`$ years and its wavelength is 3.11 kpc. For the upper point, $`\zeta _{\mathrm{node}}=30`$ ($`z_{\mathrm{node}}`$=5 kpc), the corresponding values change to $`3.4\times 10^7`$ years and 3.43 kpc, respectively. From the figure, it is clear that the maximum growth rate above $`z3`$ kpc is less sensitive to the position of the nodal point. This is because the gravitational acceleration (see eq. ) already reaches its maximum value, $`8\times 10^9`$ cm s<sup>-2</sup> at about $`3`$ kpc. Therefore, as the nodal point goes to positions higher than $`\zeta _{\mathrm{node}}=30`$, the growth time converges to $`3\times 10^7`$ years, which can be regarded as the minimum growth time of the Parker instability in the thick gaseous disk (obviously, at these heigths the gravitational force also has a non-neglegible radial component, and we are near the limit of validity of our 2D analysis). In the following section we address the structure of the final equilibrium state.

## 3 Two-dimensional Equilibria of the Undular Instability

### 3.1 Magnetohydrostatic Equations

The MHS equations are obtained by setting $`𝐯=0`$ and dropping the time-derivative terms in the MHD equations (1-3). Hence, the continuity and induction equations are of no use in this case. Due to this reason, the number of unknowns ($`\rho `$, $`B_y`$, and $`B_z`$) is larger than the number of equations (the $`y`$ and $`z`$ components of the momentum equation), and one requires an additional expression. Closure is granted with flux freezing, which results in conservation of the mass-to-flux ratio in a flux tube.

The details for the derivation of the final equilibrium states are given in Mouschovias (1974), and are summarized in Spitzer (1978). Given the magnetic vector potential $`𝐀=\widehat{e}_xA(y,z)`$

$$𝐁=\times 𝐀,$$
(21)   
and the gravitational potential

$$\psi =_0^zg(z)𝑑z,$$
(22)   
the final magnetic equilibrium is given by

$$^2A=4\pi \frac{dq}{dA}\mathrm{exp}\left(\frac{\psi }{a^2}\right).$$
(23)   
The function $`q\rho a^2\mathrm{exp}(\psi /a^2)`$ is a constant along a line of force and is given by

$$q(A)=\frac{a^2}{2}\frac{dm}{dA}\left\{_0^{\lambda _y/2}𝑑y\frac{z(y,A)}{A}\mathrm{exp}\left[\frac{\psi (y,A)}{a^2}\right]\right\}^1,$$
(24)   
where $`\lambda _y`$ is the perturbation wavelength along the initial magnetic field, and $`dm/dA`$ is the mass-to-flux ratio. As stated above, for flux freezing conditions the mass between two field lines is conserved and the mass-to-flux ratio remains constant during the evolution (i.e., is a constant of motion). Then, this ratio is determined from the initial equilibrium configuration

$$\frac{dm}{dA}=\lambda _y\frac{\rho _0(A)}{B_0(A)},$$
(25)   
where $`\rho _0(A)`$ and $`B_0(A)`$ represent the initial distributions of the density and magnetic field as functions of $`A`$, respectively.

### 3.2 Final Equilibrium States

Now, after setting the mass-to-flux ratio, one can solve equations (23) and (24) simultaneously. Following the detailed procedure described in Appendix C of Mouschovias (1974), we solve these equations by iteration. In contrast to the original work of Mouschovias, who used a constant gravity with a discontinuity at midplane, we use a smooth and continuous gravity function (eq. ). The discontinuity prevents midplane gas crossings, and he found the final equilibria of the MS modes only. We do not have such a discontinuity and are able to derive the final equilibria of both the MS and MA modes.

The initial equilibrium distributions for the density and field lines are plotted in Figure 4a. Colors are mapped from red to violet as the density decreases. The white lines represent the $`𝐁`$-field lines, and they are chosen in such a way that the magnetic flux between two consecutive lines is the same. The length scales are normalized with the effective scale height, $`H_{\mathrm{eff}}`$. First, to derive the final state of the MS mode, we added a MS perturbation to the magnetic vector potential in the initial equilibrium state,

$$\delta A(y,z)=A_0(z)C\mathrm{cos}(\frac{2\pi y}{\lambda _y})\mathrm{sin}(\frac{2\pi z}{\lambda _z}),$$
(26)   
where $`C=0.01`$ is the amplitude of the perturbation, and $`\lambda _z=2z_{\mathrm{node}}`$ (the MS mode has zero amplitude at $`z=0`$ and $`z=z_{\mathrm{node}}`$, so $`z_{\mathrm{node}}`$ corresponds to a half of the wavelength value of the principal mode along the $`z`$-axis). Figure 3 shows that when the first nodal point from midplane is $`9H_{\mathrm{eff}}`$, the most unstable horizontal wavelength is $`18H_{\mathrm{eff}}`$. Thus, we use the pair of unstable wavelengths $`(\lambda _y,\lambda _z)=(18H_{\mathrm{eff}},18H_{\mathrm{eff}})`$, and get the final equilibrium state displayed in Figure 4b. The actual computational domain for this symmetric case is $`0y9H_{\mathrm{eff}}`$ and $`0z9H_{\mathrm{eff}}`$, but for visualization purposes we extend the domain eight times in the figure. The condensations formed in the magnetic valleys and voids in the arches are clearly seen in the figure. Due to the condition imposed at midplane, $`\delta A=0`$, the field line at $`z=0`$ is not deformed at all.

For the MA case, we perturb the initial state with the perturbation,

$$\delta A(y,z)=A_0(z)C\mathrm{cos}(\frac{2\pi y}{\lambda _y})\mathrm{cos}(\frac{2\pi z}{\lambda _z}),$$
(27)   
where we now use $`\lambda _z=4z_{\mathrm{node}}`$ (in contrast to the MS mode, the MA mode has maximum amplitude at $`z=0`$ and then requires twice the wavelength value along the $`z`$-axis). As stated before, the dispersion relations shown are not sensitive to the midplane boundary conditions, and the most unstable horizontal wavelength is the same for both the MS and MA modes. Using the same nodal point as before, $`z_{\mathrm{node}}=9H_{\mathrm{eff}}`$, we set the pair of unstable wavelengths to $`(\lambda _y,\lambda _z)=(18H_{\mathrm{eff}},36H_{\mathrm{eff}})`$, and the final state of the MA mode is plotted in Figure 4c. The computational domain is now $`0y9H_{\mathrm{eff}}`$ and $`9H_{\mathrm{eff}}z9H_{\mathrm{eff}}`$, and covers the upper and lower hemispheres. For a better visual impression we extend it by a factor of four in Figure 4c. As before, one can also see condensations and voids in the figure, but their positions are now alternated between the upper and lower hemispheres. Hence the distance between successive condensations is a half of the horizontal wavelength. Also, the $`𝐁`$-field line at midplane is now undulated, with locations above and below $`z=0`$, as is characteristic of the MA mode.

The density enhancements produced by the gas that has been sliding into the magnetic valleys can be obtained from the column density of the final state. At any given location $`y`$, the final column density is

$$N_f(y)=_{z(y,A_i[z=0])}^{9H_{\mathrm{eff}}}\rho _f(y,z)𝑑z,$$
(28)   
where the subscripts $`i`$ and $`f`$ denote the initial and final states, respectively. The lower limit of the integral corresponds to the final $`z`$-coordinate of the magnetic field line that was initially located at midplane, and is labeled with $`A_i(z=0)`$. Thus, the lower limit is exactly equal to zero for the MS modes, but it is different from zero for the MA modes. In Figure 5 we plot the column density distributions for both modes (normalized to the initial column density, $`N_i=_0^{9H_{\mathrm{eff}}}\rho _i(z)𝑑z`$). The figure reveals that the MA modes drive more gas into the magnetic valleys than the MS cases. This is because the MA perturbations can gather more mass in the condensations by bending the midplane.

Another interesting quantity is the ratio of the magnetic-to-gas pressures, $`\alpha =B^2/(8\pi a^2\rho )`$. The value of this ratio varies with time and $`z`$-location, and Figure 6 shows the corresponding distributions at the initial and final equilibrium states. The initial state is plotted as a solid line, and is labeled as $`\alpha _i(z)`$. The distribution for the final MS state is plotted with dashed lines at two different $`y`$-positions: the distribution at $`y=0`$, corresponding to the center of the condensation (the magnetic valley), is labeled as $`\alpha _f(0,z)`$, and the distribution at $`y=9H_{\mathrm{eff}}`$ (the central part of the magnetic arch; see Fig. 4b) is labeled as $`\alpha _f(9,z)`$. Finally, the distribution for the final MA state at the position $`y=0`$ is shown with a dotted line, and is also labeled $`\alpha _f(0,z)`$ (the lower part of the disk contains the central part of a condensation, and the upper part of the disk has the maximum of the magnetic arches; see Fig. 4c). The $`\alpha `$-distribution of the initial state shows that our disk model is mainly supported by gas pressure near the midplane, and by magnetic pressure at high latitudes (see also Fig. 2). The distribution, however, is completely modified at the final equilibrium stages. The gas pressure increases at the condensations and the resulting pressure ratios, for both the MS and MA modes, become smaller than the initial $`\alpha `$ values. In contrast, at the voids, where the magnetic energy becomes dominant, $`\alpha `$ reaches values as high as $`10^4`$. This is because the gas is efficiently drained down from the magnetic arches, as already pointed out by Mouschovias (1974) for the case of a thin gaseous disk in a uniform gravity.

As a final comment of this section, we add that the galactic system seems to prefer the lower energy state of the MA mode. This is not apparent from the dispersion relations, which are degenerate for the MS and MA modes, but it appears in detailed numerical MHD experiments performed with the same thick disk model considered in this study. The results of these numerical simulations will be reported in a separate paper (Santillán et al. 1999b), and here we only mention the relevant result. The runs are started with the initial equilibrium state described §2.2 and, as expected, the early linear phase of the experiments follows the rates and wavelengths derived in the present linear analysis for any of the two parity modes. We also performed several experiments with random velocity perturbations, without any preferred parity or wavelength. These random perturbation experiments eventually evolve into the MA configuration, similar to the one shown Figure 4c, indicating that the MA mode is preferred over the MS mode.

## 4 Summary and Discussions

Here we have presented the linear perturbation analysis of a magnetized and warm thick disk. The dispersion relations for the undular mode of the Parker instability are derived, along with the resulting final equilibrium states. The initial disk parameters are taken from the observed distributions at the solar circle, and we assume that the gas is isothermal and that the initial field lines are parallel to the disk. Given the complexities inherent to trying to model the dynamical effects of cosmic rays, its pressure is not explicitly included in here. The resulting multi-component gaseous disk model, then, has a thermal-to-magnetic pressure ratio that decreases with $`z`$-location, and is Parker unstable. The properties of the unstable modes for five different nodal points are given in Figure 3. These nodal points correspond with the assumed extension of the disk above midplane (for $`H_{\mathrm{eff}}166`$ pc, the five cases in Figure 3 represent a disk extending up to 1.5, 2, 3, 4 and 5 kpc, respectively).

The value of the critical wavelength for the instability depends on the location of the chosen nodal point and, for disks extending between 5 and 1.5 kpc above midplane, it increases from 1.5 to 1.8 kpc, respectively. The wavelength of the fastest growing mode, however, is almost independent of the assumed nodal point, and is about 3 kpc for all the cases considered. The corresponding growth time scales are slightly more sensitive to the nodal point and, for disks extending between 1.5 and 5 kpc above midplane, the time scales decrease from 6.2 to 3.4$`\times 10^7`$ year, respectively. The minimum growth time then converges to $`3\times 10^7`$ year as the nodal point tends to large $`z`$-values. Thus, for average conditions in the solar neighborhood, the linear analysis in 2D indicates that the preferred wavelength is about 3 kpc, and that the gas condensations are formed in time scales of the order of about $`3\times 10^7`$ year. The wavelength values are larger, by a factor of about 8, than those derived for the thin disk cases, but the corresponding time scales are larger by only a factor of about 2 (Kim & Hong 1998). These are substantial differences, and indicate that the multi-component structure of the disk play an important role in the large-scale stability and evolution of the ISM.

The densities for the final equilibrium stages, on the other hand, are larger for the MA modes. The resulting MA column densities at the condensations are increased by a factor of about 3 with respect to the value of initial equilibrium stage. Now, the spiral density wave can trigger the Parker instability in the model considered here (Martos & Cox 1994), and the contrast obtained is similar to the expected density contrast between arm-interarm regions for strong waves (Elmegreen 1991). For comparison, the final-to-initial column density ratio of the fastest growing mode in a thin disk model is of about 1.2 only (Mouschovias 1974). Thus, the gas from the extended gas layers participating in the instability contribute with a fraction of about 2/3 to the total mass gathered in the condensations.

The role of self-gravity and differential rotation of the Galaxy are not included in the present study. Self-gravity may not be important at the early linear stages of the instability (e.g. Hanawa, Nakamura & Nakano 1992), but it will lead to more compact and denser condensations at the non-linear phases. Galactic differential rotation, on the other hand, has an influence at several stages of the Parker instability (e.g. Shu 1974; Zweibel & Kulsrud 1975; Balbus & Hawley 1991; Foglizzo & Tagger 1994, 1995). For instance, if the radial differential force is strong enough, a transient shearing instability also appears, and the combined Parker-shearing instability could lead to angular momentum transfer and dynamo action in disks. In the case of the 2-D undular perturbations considered in this paper, the lateral motions of the flows should be affected by the Coriolis force. If we, however, include the ignored third dimension (the radial direction) in our analysis, the vertical motion of the mixed mode with a smaller wavelength along the radial direction dominates the system, and the effects of rotation are severely reduced during the linear growth. Nonetheless, as stated by the referee, the stabilizing effects of rotation may be important at the final equilibrium stages. These are important issues that require detailed three dimensional studies with differential rotation, and should be addressed in future studies.

If the assumptions of the present work are valid, the range of growth rate values are marginally consistent with those required for the formation of giant molecular clouds in our Galaxy (e.g. Blitz & Shu 1980). Also, the most unstable wavelength in our model is somewhat larger than the corrugation distance derived by Alfaro et al. (1992) for the Carina Arm (2.4 kpc), but the condensations formed by the odd parity mode of the instability may well be associated with the origin of this observed structure. A more detailed study is required to properly address this issue, and important caveats should be borne in mind regarding the applicability of the present results. One is that of the randomness of the Galactic magnetic field topology at the kpc length scale, not included in the present modeling. Another one is the largely unknown temperature structure of the halo, and the filling factors of the different gas components. Models built from the same density and gravity distributions, but in which the magnetic field distribution is prescribed from the Galactic synchrotron emission (e.g. Martos & Cox 1998), require thermal dominance at high $`[z]`$ and are therefore Parker stable.

This leads us to a final important question if the isothermal disk assumption represents a fair description of the actual gaseous disk in our Galaxy. Here we do not differentiate between the thermal and kinetic pressures, and both are gathered in a single isothermal term with sound speed similar to the velocity dispersion of the main components extending up to $`1.5`$ kpc from midplane (Boulares & Cox 1990). Such an isothermal condition, then, can be considered as a reasonable approximation for the regions located between 1 to 1.5 kpc from the midplane. The existence of a few “anomalous” velocity components within 2 kpc (e.g. Kulkarni & Fich 1985; Reynolds 1985), and gas with a large velocity dispersion at $`z`$ of about 4 kpc (Kalberla et al. 1998), already hint that the effective sound speed should be increased somewhere within 1 and 2 kpc. The details for such a variation are presently unknown, but we are currently investigating the effects of some reasonable velocity distributions. Obviously, the loss of magnetic support provides a stabilizing effect, but the present restrictions do not indicate that the instability can be completely suppressed. A detailed discussion of the range of velocity dispersion variations and the resulting unstable mode values will be presented elsewhere.

###### Acknowledgements.

It is a big pleasure to thank Emilio Alfaro, Don Cox, Gene Parker, and Dongsu Ryu for many stimulating and informative discussions during the development of this project. We are grateful to Thierry Foglizzo, the referee, and Steve Shore, the editor, for several constructive comments. JF thanks the Korea Astronomy Observatory and Seoul National University for their warm hospitality. JF, MM and AS acknowledge partial support by DGAPA-UNAM grant IN130698, CONACyT grants 400354-5-4843E and 400354-5-0639PE, and by a R&D CRAY Research grant. The work of JK was supported by the Office of the Prime Minister through Korea Astronomy Observatory grant 99-1-200-00, and he also acknowledges the warm hospitality of the Instituto de Astronomía-UNAM. The work of SSH was supported in part by a grant from the Korea Research Foundation made in the year 1997.      

Figure Captions