text:
# Three-Dimensional Simulations of Parker, Magneto-Jeans, and Swing Instabilities in Shearing Galactic Gas Disks

## 1 Introduction

Nearly half of the interstellar medium (ISM) in the inner part of the Milky Way is estimated to be in the molecular component (e.g., Dame 1993), with the largest portion of the molecular ISM in giant molecular clouds or cloud complexes (GMCs) of masses $`10^410^6\mathrm{M}_{}`$ (Dame et al., 1987; Solomon et al., 1987), where most of star formation occurs (e.g., Williams & McKee 1997). GMCs are generally turbulent, self-gravitating, and magnetized (e.g., Blitz 1993; Crutcher 1999). In external galaxies, they often appear in clusters, forming giant molecular associations (GMAs) (e.g., Grabelsky et al. 1987; Vogel et al. 1988; Rand & Kulkarni 1990; Rand 1993; Sakamoto et al. 1999), or within the boundaries of H I superclouds (e.g., Elmegreen 1995). These giant clouds in spiral galaxies are mostly associated with spiral arms, but they sometimes appear in the interarm regions (e.g., Solomon et al. 1985; Rand 1993; Kenney 1997; Heyer & Terebey 1998). By regulating the rate and mode of star formation in disk galaxies, GMCs play a fundamental role in controlling galactic evolution and may be a key to understanding the nature of the Hubble sequence. But how GMCs form is still not clearly understood.

Traditionally, proposed mechanisms for giant cloud formation fall into two categories: stochastic coagulation of smaller clouds (following Oort 1954; e.g., Kwan 1979) or collective effects involving instabilities (see e.g., Elmegreen 1996). Existing work suggests that collisional agglomeration may proceed too slowly (e.g., Scoville & Hersh 1979; Elmegreen 1990, 1995). Collisions in most cases result in disruption rather than merger (e.g., Lattanzio et al. 1985; Kim et al. 1999; Klein, Woods, & McKee 2001); the available mass in smaller clouds is not enough to form GMCs through collisional build-up (Heyer & Terebey, 1998; Blitz & Williams, 1999); collisional agglomeration would lead to magnetically supercritical clouds and active star formation before GMC-scale masses are reached (Ostriker, Gammie, & Stone, 1999; Ostriker, Stone, & Gammie, 2001); GMC chemistry may be inconsistent with an extended build-up phase (e.g., van Dishoek et al. 1993). All of these and other difficulties point to the need for collective effects (Blitz & Shu, 1980). This second means of building GMCs and GMAs involves large-scale dynamical instabilities that coherently amass material over scales large compared to the sizes of diffuse atomic clouds (see e.g., Elmegreen 1995). The small-scale structure of the medium plays a secondary role, with the modal gathering of material operating collectively on the whole distribution. Coherently-growing unstable modes that produce condensations depend on various effects, including self-gravity, the Coriolis force, sheared azimuthal motion, magnetic tension, and magnetic buoyancy, all of which have been shown to interact in unstable fashion.111Thermal instabilities involving changes in the microphysical kinetic temperature are important in forming cold clouds on small scales (Field 1965; see also e.g., Hennebelle & Pérault 1999 and Burkert & Lin 2000 for recent simulations); if dissipation of kinetic energy in collisions of cloudlets acts effectively like microscopic cooling for the turbulent “cloud fluid”, then analogous mechanisms may also be important on large scales (e.g., Elmegreen 1989, 1991).

The two most well-known condensation instabilities are the Parker instability (Parker, 1966), in which buoyancy causes dominantly in-plane magnetic fields to buckle, with matter collecting in valleys; and the swing amplifier, in which epicyclic motion conspires with galactic shear to enable self-gravitating enhancement of overdense regions (Goldreich & Lynden-Bell, 1965b; Toomre, 1981). Another condensation instability that has been shown to be efficient in low-shear, magnetized regions is the magneto-Jeans instability (MJI), in which magnetic tension of in-plane fields counters the stabilizing effect of the Coriolis force, allowing self-gravitating contraction of overdense regions (Elmegreen, 1987; Kim & Ostriker, 2001). Because the growth rates of Parker and MJI modes increase with magnetic field strength and density, the compression of the ISM induced by passage through stellar spiral arms may be important in triggering GMC and GMA formation (e.g., Mouschovias, Shu, & Woodward 1974; Blitz & Shu 1980). Observed distributions of H II regions and OB star complexes in a “beads on a string” pattern along spiral arms may be the direct consequence of cooperation among the Parker instability, Jeans instability, and spiral arm potentials (e.g., Elmegreen & Elmegreen 1983).

In Kim & Ostriker (2001, hereafter Paper I), we presented studies of detailed nonlinear evolution of self-gravitating modes in shearing, razor-thin galactic disks. We showed that formation of bound condensations via the swing instability is subject to “threshold” behavior, with models having sufficiently small values of the Toomre $`Q`$ parameter (see eq.  below) producing bound clouds, and models with larger values of $`Q`$ remaining nonlinearly stable. We found that the threshold value of $`Q`$ is in the range $`1.21.4`$ (for magnetic field strengths from zero up to thermal-equipartition values), consistent with an apparent mean threshold $`Q_{\mathrm{th}}=1.4`$ for star formation evident from observational studies (Kennicutt, 1989; Martin & Kennicutt, 2001). In Paper I, we also elucidated the role of MJI in forming bound condensations by showing that in low-shear regions where the swing amplifier cannot operate, the mechanism that leads to growth of perturbations is quite different from that in high-shear regions, and that the presence of magnetic fields is crucial to the instability. In the absence of magnetic fields, specific vorticity is conserved; contraction of vorticity-conserving masses leads to increased Coriolis forces that exceed the increase in self-gravity such that further contraction is prevented (in the absence of a shearing background). In magnetized systems, however, tension forces can transfer vorticity from one region to another, such that condensation driven by self-gravity can proceed.

Effects of spiral arm gravity on the growth of giant clouds were explicitly addressed by Kim & Ostriker (2002, hereafter Paper II) using two-dimensional models. Paper II showed that for a spiral arm inside corotation, the enhanced surface density and shear gradients in the arm region promote the growth of gaseous spurs. These spurs jut nearly perpendicularly from the outer (downstream) side of the arm, and trail in the same sense as the arm itself; this is consistent with observed spur-like dust lanes in M51 (Scoville et al., 2001). We argued that the formation of spurs can be understood as a natural consequence of the MJI within spiral arms, where compression leads to a reduced (or negative) shear gradient that allows these self-gravitating modes to grow over an extended time. Paper II showed that bound clouds of typical mass $`4\times 10^6\mathrm{M}_{}`$ can form from spur fragmentation; such condensations potentially represent structures that would initially be observed as GMAs, and subsequently evolve to become bright OB associations.

Because of their restriction to razor-thin geometry, our previous two-dimensional models could not capture the potential consequences of the Parker instability on cloud formation. While previous work on the Parker instability has studied the weak-shear limit (see linear-theory analysis of Shu 1974, Foglizzo & Tagger 1994, 1995) and the case of solid-body rotation (see linear theory of Zweibel & Kulsrud 1975 and simulations by Chou et al. 2000 and Kim, Ryu, & Jones 2001), and indeed shown that rotation and shear tend to be stabilizing, no previous study focused primarily on Parker modes has incorporated shear at a realistic “average” galactic level. Although the Parker instability has been proposed to prompt the formation of molecular clouds (e.g., Mouschovias, Shu, & Woodward 1974; Blitz & Shu 1980), existing numerical simulations indicate that the Parker instability alone is unable to produce structures of high enough surface density to represent GMCs (Basu, Mouschovias, & Paleologou, 1997; J. Kim et al., 1998, 2000; Santillán et al., 2000; Kim, Ryu, & Jones, 2001).222In addition, a random component of galactic magnetic fields whose strength is almost comparable to the mean uniform component is shown to suppress the Parker instability significantly (Kim & Ryu, 2001).

To study the coupling of the Parker instability with self-gravitating modes, Elmegreen (1982a, b) performed linear stability analyses for combined Parker-Jeans modes and concluded that clouds of mass $`10^510^6\mathrm{M}_{}`$ form in the density regime for which self-gravity significantly modifies the Parker instability; this regime of densities is expected to be found in spiral arm regions only. Hanawa, Nakamura, & Nakano (1992) included solid-body rotation in a linear-theory analysis of the Parker-Jeans instability, and concluded that the transfer of angular momentum out of growing condensations may be a limiting factor. This transfer of angular momentum may be accomplished by magnetic stresses much as it occurs in the two-dimensional MJI. A nonlinear simulation of the Parker-Jeans instability in a disk with solid-body rotation was among the models performed by Chou et al. (2000), who found that filamentary structures tend to form perpendicular to the mean magnetic field.

To our knowledge, there have previously been neither linear-theory nor numerical studies that incorporate all of the aforementioned elements that can contribute to three-dimensional instabilities in galactic disks: buoyant (horizontal) magnetic fields, self-gravity, rotation, and shear. In this paper, we investigate nonlinear evolution of galactic disks allowing for all these fundamental physical agents. This work extends our previous studies (Paper I) of structure formation via gravitational instability by including the effects of the Parker instability. Our primary objectives are to study how the Parker instability develops in the presence of strong shear, and to compare the eventual evolution of disks under the various instabilities separately or in combination. We can thus assess the potential for gravitationally-bound cloud formation from a variety of dynamical mechanisms. Since vertically-integrated, razor-thin disks are known to overestimate self-gravity at the midplane (e.g., Toomre 1964), we will also study analytically and numerically how the stability of three-dimensional disks to in-plane self-gravitating modes changes as their thickness varies.

In our model, disks are local (see §2.1), isothermal in both space and time, and vertically stratified. Initial magnetic fields are plane-parallel, pointing in the azimuthal direction. To model realistic conditions in disk galaxies, we include both external gravity arising from stars (as a time-independent potential) and gaseous self-gravity, which together counterbalance thermal and magnetic pressure gradients in an unperturbed state. We do not, however, consider in this paper other features such as stellar spiral arms and effects of cosmic rays that may enhance the magnetic buoyancy force. We begin by solving the governing equations to obtain initial static equilibrium density distributions, and then explore their gravitational stability to axisymmetric perturbations. We next select several parameter sets that represent “average disk” or “spiral arm” regions. After applying small-amplitude perturbations to the corresponding density profiles, we follow the nonlinear dynamical development with three-dimensional direct numerical simulations.

The paper is organized as follows: In §2, we describe the basic MHD equations in the local “shearing-sheet” approximation, the computational methods we use, and our model parameters. In §3, we present initial magnetohydrostatic equilibria in the presence of both self-gravity and external gravity, and analytically examine their axisymmetric stability by making simplifying assumptions. In §4, we present the results of three-dimensional simulations for the Parker instability without self-gravity (§4.1), generalized (Parker/swing) gravitational instability in shearing regions (§4.2), and MJI in nonshearing regions (§4.3). We summarize our results and discuss the implications of present work for galactic cloud formation in §5.

## 2 Numerical Methods and Model Parameters

### 2.1 Basic Equations and Numerical Methods

In this paper, we study nonlinear evolution of vertically stratified, differentially rotating, self-gravitating galactic gas disks with threaded magnetic fields in the azimuthal direction. We solve fully three-dimensional, compressible, ideal MHD equations in a local Cartesian reference frame whose center lies at a galactocentric radius $`R_0`$ and orbits the galaxy with a fixed angular velocity $`\mathrm{\Omega }\mathrm{\Omega }(R_0)`$. In this local frame, radial, azimuthal, and vertical coordinates are represented by $`xRR_0`$, $`yR_0(\varphi \mathrm{\Omega }t)`$, and $`z`$, respectively, and terms associated with the curvature effect are neglected (Goldreich & Lynden-Bell, 1965b; Julian & Toomre, 1966). The local-frame equilibrium background velocity relative to the center of the box at $`x=y=z=0`$ is then reduced to $`𝐯_0q\mathrm{\Omega }x\widehat{𝐲}`$, where

$$q\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}R}|_{R_0},$$
(1)
is the local dimensionless shear rate. In terms of $`q`$, the local epicyclic frequency $`\kappa `$ is given by $`\kappa ^2(42q)\mathrm{\Omega }^2`$. The resulting “shearing-sheet” equations expanded in the local frame are

$$\frac{\rho }{t}+(\rho 𝐯)=0,$$
(2)
$$\frac{𝐯}{t}+𝐯𝐯=\frac{1}{\rho }P+\frac{1}{4\pi \rho }(\times 𝐁)\times 𝐁+2q\mathrm{\Omega }^2x\widehat{𝐱}2𝛀\times 𝐯\mathrm{\Phi }_s+𝐠_{\mathrm{ext}},$$
(3)
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁),$$
(4)
$$^2\mathrm{\Phi }_s=4\pi G\rho ,$$
(5)
and

$$P=c_s^2\rho ,$$
(6)
(cf, Hawley, Gammie, & Balbus 1995; Stone et al. 1996). Here, $`𝐯`$ is the velocity in the local frame, $`\rho `$ is the mass density, $`P`$ is the thermal pressure, $`𝐁`$ is the magnetic field, $`\mathrm{\Phi }_s`$ is the self-gravitational potential, $`𝐠_{\mathrm{ext}}`$ is the external gravitational acceleration, $`G`$ is the gravitational constant, and $`c_s`$ is the isothermal sound speed. In equation (3), the third and fourth terms on the right hand side represent the tidal and Coriolis forces, respectively, which balance each other in the initial equilibrium with $`𝐯=𝐯_0`$. As expressed by equation (6), for simplicity we adopt an isothermal equation of state for all the simulations presented in this paper.

An initial vertical equilibrium is defined by the balance between the total (thermal plus magnetic) pressure gradient and the total (self plus external) gravitational force. For $`𝐠_{\mathrm{ext}}`$ representing the external gravity from stars, we adopt the simple form

$$𝐠_{\mathrm{ext}}=\frac{g_{}}{H_{}}𝐳,$$
(7)
where $`H_{}`$ is the scale height of stars and $`g_{}`$ is a free parameter. Since most of gas is distributed within $`|z|<H_{}`$ and the small-$`|z|`$ part of the realistic gravity deduced from stellar tracer populations can be described by a linear function of height (e.g., Bahcall 1986; Kuijken & Gilmore 1989), $`𝐠_{\mathrm{ext}}`$ given in equation (7) is a fairly good approximation in studying dynamics of the gaseous component. Although a linear gravity law is known to strengthen the Parker instability slightly compared to the case with a more realistic gravity law, the growth rate of the most unstable mode differs by only a factor of 1.2 (Kim & Hong, 1998). As we shall show in §3, the total gravity that drives the Parker instability in our models is similar to the external gravity adopted by other authors (Giz & Shu, 1993; Kim & Hong, 1998).

The numerical solutions to the dynamical equations (2)$``$(7) are obtained using a modified version of the ZEUS code originally developed by Stone & Norman (1992a, b). ZEUS is a time-explicit, operator-split, finite-difference method for solving the MHD equations on a staggered mesh, and has been demonstrated to provide robust solutions to a wide variety of astrophysical problems. ZEUS employs “constrained transport” to guarantee that $`𝐁=0`$ within machine precision, and the “method of characteristics” for accurate propagation of Alfvénic disturbances (Evans & Hawley, 1988; Hawley & Stone, 1995). For less diffusive advection of hydrodynamic variables, we apply a velocity decomposition method which treats the contribution from the background shearing parts as source terms, while solving the transport step using only the perturbed parts. In order to minimize the errors arising from discontinuities in the flow characteristics across the $`x`$-boundaries, the ghost zones adjoining the $`x`$-boundaries are kept active (Paper I). In the $`\widehat{x}`$ and $`\widehat{y}`$ directions, we implement shearing box boundary conditions in which the $`x`$-boundaries are shearing-periodic and the $`y`$-boundaries are perfectly periodic (Hawley, Gammie, & Balbus, 1995). In the $`\widehat{z}`$ direction, we adopt free (outflow) boundary conditions which set all the variables in the boundary zones equal to the corresponding values in the adjacent active zones (Stone & Norman, 1992a). These boundary conditions permit outgoing waves and thus physical quantities such as mass, momentum, and magnetic flux can flow out of the simulation domain. Since density is very low near the vertical boundaries, however, simulation results generally show that the total mass is conserved within 1%.

For computation of the self-gravitational potential, we employ a hybrid technique that combines the Green function method to compute the gravitational kernels (Miyama, Narita, & Hayashi, 1987) and the Fourier transform method in shearing-sheet coordinates (Gammie, 2001). Our high-resolution models have a $`256\times 128\times 128`$ grid for ($`x`$, $`y`$, $`z`$); we find that the Poisson solver typically takes about 48% of the total cpu time. We have checked our implementation of the ZEUS code on various test problems including two-dimensional gravitational instability and propagation of gravito-MHD waves in three dimensions.

### 2.2 Model Parameters

Our model disks are highly-idealized, smoothed versions of what is physically a much more complex, multiphase interstellar medium (e.g., Field, Goldsmith, & Habing 1969; Cox & Smith 1974; McKee & Ostriker 1977; Heiles 2001). We adopt an isothermal condition in both space and time, with an effective isothermal speed of sound $`c_s=7\mathrm{km}\mathrm{s}^1`$, corresponding to a mean Galactic thermal pressure $`P/k20004000\mathrm{K}\mathrm{cm}^3`$ (Heiles, 2001) and mean midplane density $`n_\mathrm{H}0.6\mathrm{cm}^3`$ (Dickey & Lockman, 1990); external galaxies have similar gaseous velocity dispersions. One of the parameters in our simulations is the rotation rate $`\mathrm{\Omega }`$. When reporting dimensional results, we normalize to the solar neighborhood value $`\mathrm{\Omega }=26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1`$ (Binney & Tremaine, 1987). For a flat rotation curve, $`q1`$ and thus $`\kappa 2^{1/2}\mathrm{\Omega }`$; we allow for varying $`q`$ to represent both “average” regions and reduced-shear (spiral arm) regions.

Our initial disk models are vertically stratified, with initial density distribution $`\rho _0(z)`$. The corresponding surface density is given by $`\mathrm{\Sigma }_0_{\mathrm{}}^{\mathrm{}}\rho _0(z)𝑑z`$. To characterize $`\mathrm{\Sigma }_0`$ in a dimensionless fashion, we use the Toomre $`Q`$ stability parameter

$$Q\frac{\kappa c_s}{\pi G\mathrm{\Sigma }_0}.$$
(8)
For an infinitesimally-thin disk, the critical value for axisymmetric gravitational stability is $`Q_c=1`$ (Toomre 1964; see also e.g., Binney & Tremaine 1987; Shu 1992). We shall show in §3 that the finite thickness of disks reduces $`Q_c`$ (cf. Goldreich & Lynden-Bell 1965a). For local gas surface density of $`\mathrm{\Sigma }_0=13\mathrm{M}_{}\mathrm{pc}^2`$ comparable to the solar neighborhood (cf, Paper I; Kuijken & Gilmore 1989), $`Q1.5`$; we adopt this value of $`Q`$ for our “average disk” fiducial value.

The initial magnetic field in our models points in the azimuthal ($`\widehat{𝐲}`$) direction and varies with height as $`𝐁_0=B_0(z)\widehat{𝐲}`$; its magnitude is parameterized by

$$\beta \frac{c_s^2}{v_\mathrm{A}^2},$$
(9)
where the Alfvén speed is $`v_\mathrm{A}B_0/\sqrt{4\pi \rho _0}`$. As an initial condition, we fix $`\beta =1`$ spatially for all the MHD models in this paper, so that initially the Alfvén speed is the same as the isothermal sound speed everywhere. This implies $`B_0(z)\sqrt{\rho _0(z)}`$.

We explore the nonlinear evolution of disks under both average disk and other galactic conditions by changing $`\mathrm{\Sigma }_0`$ and $`q`$. In order to model a portion of a gas disk inside a spiral arm, we use $`\mathrm{\Sigma }_0=65\mathrm{M}_{}\mathrm{pc}^2`$, which is five times larger than the average disk value. Since $`\kappa \mathrm{\Sigma }^{1/2}`$ and $`Q\mathrm{\Sigma }^{1/2}`$ for isothermal gas from the constraint of potential vorticity conservation (Balbus & Cowie 1985; Paper II), this higher surface density corresponds to $`Q0.7`$, which we adopt as a fiducial “spiral arm” value. We allow the shear parameter $`q`$ to vary from model to model, but fix it spatially in a given model. Note that this arbitrary implementation of background shear is not consistent with realistic shear that varies rapidly across the spiral arm (Paper II). Nevertheless, our high surface density models allow us to focus on the effects of enhanced self-gravity inside the spiral arm without the additional complications of in-plane variations of the background.

We parameterize the magnitude of the external gravity by

$$s_0\frac{4\pi G\rho _{\mathrm{}}H_{}}{g_{}},$$
(10)
which measures the relative strength of self-(gaseous) gravity and external (stellar) gravity in the equilibrium state (cf, Elmegreen 1982a). Here, $`\rho _{\mathrm{}}`$ denotes the midplane density in the absence of both magnetic field and external gravity (i.e., $`\beta =\mathrm{}`$ and $`s_0=\mathrm{}`$; see §3). If the stellar distribution approximates a vertical self-gravitating equilibrium with vertical velocity dispersion $`\sigma _{,z}`$, one can rewrite equation (10) as

$$s_0=\left(\frac{\sigma _{,z}\mathrm{\Sigma }_0}{c_s\mathrm{\Sigma }_{}}\right)^2,$$
(11)
where $`\mathrm{\Sigma }_{}`$ denotes the stellar surface density.333For an application to a central-mass dominated Keplerian disk, $`g_{\mathrm{ext}}=\mathrm{\Omega }^2z`$. Then, $`s_04\pi G\rho _{\mathrm{}}/\mathrm{\Omega }^2=2Q^2`$. Notice that $`s_0`$ is now expressed in terms of observable quantities. Since $`\mathrm{\Sigma }_{}35\mathrm{M}_{}\mathrm{pc}^2`$ and $`\sigma _{,z}20\mathrm{km}\mathrm{s}^1`$ (e.g., Kuijken & Gilmore 1989; Holmberg & Flynn 2000), $`s_01`$ is suitable for “average disk” conditions, while $`s_025`$ would describe a spiral arm region (taking the stellar spiral perturbation as weak).

Our simulation domain is a rectangular parallelepiped of size $`L_x\times L_y\times L_z`$, with a vertical range of $`|z|L_z/2=4H`$, where $`H`$ is the gaseous scale height (see eq.  below). The horizontal dimensions of the box are $`L_x=L_y=25H`$ or $`17H`$ for $`s_0=1`$ or $`s_0=25`$ models, respectively. This size corresponds to the most unstable wavelength of the axisymmetric gravitational instability in each model with $`\beta =1`$ (see §3.2). Conversion of $`H`$ into a physical value requires solving equations governing magnetohydrostatic equilibria, which will be presented in the next section. Using $`c_s/\mathrm{\Omega }=270\mathrm{pc}(c_s/7\mathrm{km}\mathrm{s}^1)(\mathrm{\Omega }/26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1)^1`$ and borrowing the results from Figure 1, we obtain for high surface-density (spiral arm) regions ($`s_0=25`$), $`H50`$ pc for $`\beta =\mathrm{}`$ and $`H80`$ pc for $`\beta =1`$ models. Under average disk conditions ($`s_0=1`$), $`H160`$ pc for $`\beta =\mathrm{}`$ and $`H210`$ pc for $`\beta =1`$ cases. The former value for an unmagnetized model is consistent with the usual inner-galaxy scale height of $``$130-170 pc derived from the hydrostatic force balance between external gravity and turbulent pressure in the average disk condition (see e.g., Boulares & Cox 1990; Lockman & Gehman 1991; J. Kim et al. 2000). On the other hand, the latter value for a $`\beta =1`$ model is similar to the direct estimate $``$220 pc for the H I scale height obtained by tangent point observations near the solar circle (Malhotra, 1995), beyond which the H I disk is known to flare dramatically (Kulkarni, Blitz, & Heiles, 1982).

The larger initial scale height of magnetized models reflect the potential importance of magnetic pressure in supporting a galactic disk against gravity. If magnetic support is significant but is neglected in observationally-based determinations of the scale height from surface density and velocity dispersion, then $`H`$ could be underestimated. We note, however, that the dynamical rearrangement of magnetic flux in our simulations always results in a reduction of the scale height over time, so that the initial values of $`H`$ in our $`\beta =1`$ models is an upper limit for a real disk of high magnetization.

The total mass $`M_{\mathrm{tot}}\mathrm{\Sigma }_0L_xL_y`$ contained in the simulation box corresponds to $`3.6\times 10^8\mathrm{M}_{}(H/210\mathrm{p}\mathrm{c})^3`$ for $`s_0=1,\beta =1`$ models, $`1.2\times 10^8\mathrm{M}_{}(H/80\mathrm{p}\mathrm{c})^3`$ for $`s_0=25,\beta =1`$ models, and $`0.5\times 10^8\mathrm{M}_{}(H/50\mathrm{p}\mathrm{c})^3`$ for $`s_0=25,\beta =\mathrm{}`$ models.

For all the simulations reported in this paper, we apply initial perturbations only to the density. We first generate a Gaussian random variable $`\delta `$ using the usual Box-Muller algorithm (cf, Press et al. 1992) and then construct its power spectrum in the Fourier space such that $`|\delta _k|^2k^{11/3}`$ for $`1kL_x/2\pi 32`$ and zero for $`32<kL_x/2\pi `$. We normalize by fixing the standard deviation of $`\delta `$ in real space to be 1%. The final form of density perturbations is given by the product of $`\delta `$ and the background density profile $`\rho _0`$.

Finally, we remark on a few dynamical timescales of note. The rotation time is defined by $`t_{\mathrm{rot}}1/\mathrm{\Omega }=3.8\times 10^7\mathrm{yrs}(\mathrm{\Omega }/26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1)^1`$. The shearing time is $`t_{\mathrm{sh}}1/q\mathrm{\Omega }`$, which is equal to $`t_{\mathrm{rot}}`$ for a flat rotation curve ($`q=1`$). Related to $`t_{\mathrm{sh}}`$ is the swing amplification time scale, which is approximately $`t_{\mathrm{sw}}(45)t_{\mathrm{sh}}`$ (cf, Toomre 1981; Paper I). The orbital period is $`t_{\mathrm{orb}}2\pi /\mathrm{\Omega }=2.4\times 10^8\mathrm{yrs}(\mathrm{\Omega }/26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1)^1`$, which is the time unit we adopt for our presentation. The typical growth time of the Parker instability is given by the Alfvén crossing time over the scale height

$$t_{\mathrm{Parker}}\frac{H}{v_\mathrm{A}}=3.0\times 10^7\mathrm{yrs}\left(\frac{H}{210\mathrm{pc}}\right)\left(\frac{v_\mathrm{A}}{7\mathrm{km}\mathrm{s}^1}\right)^1,$$
(12)
while a characteristic gravitational contraction time is

$$t_{\mathrm{grav}}\frac{c_s}{G\mathrm{\Sigma }_0}=1.2\times 10^8\mathrm{yrs}\left(\frac{c_s}{7\mathrm{km}\mathrm{s}^1}\right)\left(\frac{\mathrm{\Sigma }_0}{13\mathrm{M}_{}\mathrm{pc}^2}\right)^1,$$
(13)
(Paper I). Note that the Parker instability involves vertical motions of gas, while fluid motions in the gravitational instability are mostly parallel to the galactic plane at least in the linear regime.

## 3 Initial Equilibria and Axisymmetric Stability

### 3.1 Initial Magnetohydrostatic Equilibria

In the absence of perturbations, our isothermal disk is in vertical magnetohydrostatic equilibrium with a density profile $`\rho _0(z)`$. Thermal and magnetic pressure forces balance against the combination of gaseous self-gravity and external gravity. Stronger gravity leads to a thinner disk. In this section, we study the change of disk thickness as $`s_0`$ varies.

For simplicity, let us define $`\rho _{00}`$ be the midplane value of $`\rho _0`$ and

$$H\frac{1}{\rho _{00}}_0^{\mathrm{}}\rho _0(z)𝑑z,$$
(14)
be the gaseous scale height for general cases, while $`\rho _{\mathrm{}}`$ and $`H_{\mathrm{}}`$ are the counterparts for purely self-gravitating, unmagnetized disks (i.e., $`s_0=\mathrm{}`$ and $`\beta =\mathrm{}`$). From the definitions, it follows for fixed $`\mathrm{\Sigma }_0`$ that $`H=\mathrm{\Sigma }_0/(2\rho _{00})=H_{\mathrm{}}\rho _{\mathrm{}}/\rho _{00}`$ and using the Spitzer (1942) isothermal disk solution

$$H_{\mathrm{}}\frac{c_s}{(2\pi G\rho _{\mathrm{}})^{1/2}}=\frac{c_s^2}{\pi G\mathrm{\Sigma }_0}=280\mathrm{pc}\left(\frac{c_s}{7\mathrm{km}\mathrm{s}^1}\right)^2\left(\frac{\mathrm{\Sigma }_0}{13\mathrm{M}_{}\mathrm{pc}^2}\right)^1.$$
(15)
Arranging equations (3) and (5) – (7) for a static equilibrium in the vertical direction, we are left with

$$\frac{1}{2}\left(1+\frac{1}{2\beta }\right)\frac{d^2\mathrm{ln}\stackrel{~}{\rho }}{d\zeta ^2}=\frac{H_{\mathrm{}}}{H}\stackrel{~}{\rho }\frac{1}{s_0},$$
(16)
where the dimensionless height $`\zeta z/H_{\mathrm{}}`$ and the dimensionless density $`\stackrel{~}{\rho }(z)\rho _0(z)/\rho _{00}`$. Equation (16) is to be integrated subject to the constraint (14).

In the two limiting cases of $`s_0`$, equations (14) and (16) yield analytic solutions. For $`s_01`$ corresponding to a strongly self-gravitating disk and negligible external gravity, we obtain $`\stackrel{~}{\rho }=\mathrm{sech}^2(z/H)`$ along with $`H/H_{\mathrm{}}=1+1/2\beta `$ (e.g., Spitzer 1942; Elmegreen & Elmegreen 1978), while for $`s_01`$ corresponding to weak self-gravity compared to external gravity, $`\stackrel{~}{\rho }=\mathrm{exp}(\pi z^2/4H^2)`$ along with $`H/H_{\mathrm{}}=(\pi s_0)^{1/2}(1+1/2\beta )^{1/2}/2`$. Note that the latter Gaussian function is the usual solution for density in an equilibrium under linear gravity, with $`4H^2/\pi =2c_s^2H_{}/g_{}`$.

We solve equations (14) and (16) iteratively: for general $`s_0`$, we arbitrarily choose an initial $`H/H_{\mathrm{}}`$ and numerically integrate equation (16) to find $`\stackrel{~}{\rho }(\zeta )`$. If the calculated $`\stackrel{~}{\rho }(\zeta )`$ does not satisfy the constraint (14), we modify $`H/H_{\mathrm{}}`$ and iterate until convergence is attained. The resulting $`H/H_{\mathrm{}}`$ as a function of $`s_0`$ is shown in Figure 1 where the solid line corresponds to the unmagnetized case ($`\beta =\mathrm{}`$), while dotted and dashed lines are for sub-equipartition ($`\beta =10`$) and equipartition ($`\beta =1`$) magnetic field cases, respectively. To compare $`H`$ with a length scale independent of the gas disk’s mass, Figure 1 also plots $`H\mathrm{\Omega }/c_s=(H/H_{\mathrm{}})(Q\mathrm{\Omega }/\kappa )`$ as a function of $`s_0`$. One may think of increasing $`s_0`$ as either decreasing the stellar content ($`\mathrm{\Sigma }_{}`$) while holding the gas content ($`\mathrm{\Sigma }_0`$) fixed, or increasing $`\mathrm{\Sigma }_0`$ for a given value of $`\mathrm{\Sigma }_{}`$. Since $`H_{\mathrm{}}`$ is dependent on $`\mathrm{\Sigma }_0`$, the first meaning of $`s_0`$ is helpful in interpreting the $`H/H_{\mathrm{}}`$ curves: moving toward larger (smaller) $`s_0`$, a gas disk thickens (thins) as the stellar content relative to gas decreases (increases). The magnetic pressure support has a larger impact at higher $`s_0`$. On the other hand, since $`\mathrm{\Omega }/c_s`$ is independent of $`s_0`$ and assumed to be constant at a given radius in a galaxy, the second meaning of $`s_0`$ helps interpret the $`H\mathrm{\Omega }/c_s`$ curves: as gas moves from interarm regions into spiral arms where density and $`s_0`$ are higher, the gas layer becomes thinner in response to its own enhanced self-gravity. For instance, $`H(s_0=1)/H(s_0=25)3`$ when $`\beta =\mathrm{}`$. Note that to take into account the effects of compressed magnetic fields, $`\beta `$ should change correspondingly in spiral arms.

As mentioned before, the effective driving force, $`g_{\mathrm{eff}}`$, of the Parker instability in our models is the sum of the external gravity and the self-gravity of a disk,

$$g_{\mathrm{eff}}g_{\mathrm{ext}}\frac{d\mathrm{\Phi }_s}{dz}.$$
(17)
In Figure 2 we plot $`g_{\mathrm{eff}}`$ as solid lines together with contributions of gaseous self-gravity (dotted lines) and external stellar gravity (short-dashed lines). For average disk conditions ($`Q1.5`$ and $`s_01`$), for $`|z|/H\stackrel{<}{}1`$, external gravity is as important as self-gravity in establishing an equilibrium (Fig. 2$`a`$). For comparison, we also plot the external gravity models adopted by Giz & Shu (1993) and Kim & Hong (1998), which are in good agreement with our total average disk gravity model. On the other hand, strongly-compressed spiral arm regions are entirely dominated by self-gravity (Fig. 2$`b`$).

### 3.2 Axisymmetric Stability

Dynamics of galactic gaseous disks have often been studied by using height-integrated equations, accurate when the vertical height of the disk is much smaller than the physical scale of interest (e.g, Binney & Tremaine 1987; Shu 1992). One of the drawbacks of this thin-disk approximation is that it overestimates self-gravity for horizontal scales approaching the disk thickness (e.g., Toomre 1964). As mentioned earlier, the thin-disk approximation leads to $`Q_c=1`$ for axisymmetric gravitational instability. Goldreich & Lynden-Bell (1965a) studied the stability of uniformly rotating disks of finite thickness, and showed that $`Q_c=0.676`$ for a purely self-gravitating ($`s_0=\mathrm{}`$), isothermal disk (see also Gammie 2001). In this section, we explore the axisymmetric stability of isothermal disks under both self-gravity and external gravity.

For analytic simplicity, we make an approximation that fluid motions are all lateral, confined to the $`x`$$`y`$ plane. This ignores the small vertical oscillations of the disk as well as the potential effects of the Parker instability (which would be nonaxisymmetric anyway). We perform a standard linear stability analysis for axisymmetric modes, assuming perturbations are a linear superposition of Fourier modes, $`\mathrm{exp}(i\omega tikx)`$, of frequency $`\omega `$ and wavenumber $`k`$. For the gravitational force due to a perturbed density distribution $`\rho _1=\rho _1(z)e^{ikx}`$ with $`k0`$, the Poisson equation (5) yields an approximate solution $`d\mathrm{\Phi }_1/dx=4\pi i\mathrm{sgn}(k)G\rho _1(0)He^{ikx}/(1+|k|H)`$ at the disk midplane. This expression for the perturbed force is exact for an exponential density distribution $`\rho _1(z)e^{z/H}`$ (e.g., Elmegreen 1987), and asymptotically approaches the exact solutions for arbitrary density profiles when $`|k|H0`$ or $`\mathrm{}`$. For the vertical density profiles that we found in the previous section, the largest fractional difference between the true and approximate values for the perturbed radial gravity is only $`15\%`$, which occurs at $`|k|H1`$. The procedure for driving the dispersion relation is quite straightforward and we write down only the result,

$$\omega ^2=\kappa ^2+(c_s^2+v_\mathrm{A}^2)k^2\left(\frac{|k|H}{1+|k|H}\right)4\pi G\rho _{00}.$$
(18)
For short-wavelength perturbations $`kH1`$, equation (18) recovers the dispersion relation of waves in a rotating medium with uniform density $`\rho _{00}`$ in three dimensions (Chandrasekhar, 1961). For long-wavelength perturbations $`kH1`$, equation (18) reduces to the dispersion relation in a razor-thin disk with surface density $`\mathrm{\Sigma }_0=2\rho _{00}H`$.

Using the $`Hs_0`$ relationship that we found in the previous subsection, we compute $`k_{\mathrm{max}}`$ which minimizes the right hand side of equation (18) for given $`s_0`$ and $`\mathrm{\Sigma }_0`$ (and thus maximizes the growth rate). The results are shown in Figure 3$`a`$. The most unstable wavelength for axisymmetric perturbations is typically $`\lambda _{\mathrm{max}}=2\pi /k_{\mathrm{max}}20H`$ for average disk conditions (i.e., $`s_0=1`$) and decreases slightly as gaseous self-gravity becomes more important (larger $`s_0`$). For fixed gas content, the $`k_{\mathrm{max}}H_{\mathrm{}}`$ curves show that relatively weaker external gravity (larger $`s_0`$) requires a larger spatial scale for perturbations to become unstable. Note that $`2/H_{\mathrm{}}=2\pi G\mathrm{\Sigma }_0/c_s^2`$ is equal to the Jeans wavenumber of a razor-thin disk (cf, Paper I).

We calculate $`Q_c=\kappa (2\pi G\rho _{\mathrm{}})^{1/2}`$ corresponding to marginal stability modes ($`\omega =0`$) from equation (18) and §3.1, and plot the results in Figure 3$`b`$. Our simplified analytic approach gives $`Q_c=0.647`$ for unmagnetized disks with no external gravity ($`s_0\mathrm{}`$), which is close to $`Q_c=0.676`$ of Goldreich & Lynden-Bell (1965a)’s rigorous analysis. Critical values $`Q_c`$ for other $`s_0`$ and $`\beta `$ are shown in Figure 3$`b`$. We also perform axisymmetric ($`x`$$`z`$) simulations using our numerical code to find $`Q_c`$ for realistic density distributions obtained in §3.1. The simulation results are displayed in Figure 3$`b`$ as filled circles ($`\beta =\mathrm{}`$), triangles ($`\beta =10`$), and squares ($`\beta =1`$), which are in excellent agreement with the results of the analytic approach. For average disk conditions ($`s_0=1`$), $`Q_c`$ is found to be $``$0.75, 0.72, and 0.57 for $`\beta =\mathrm{}`$, 10, and 1 cases, respectively. These represent approximately a 25% reduction relative to the corresponding $`Q_c`$ values of razor-thin disks. This suggests that the average disk with $`Q1.5`$ is quite stable to axisymmetric gravitational instability.

We do not attempt in this paper to pinpoint threshold $`Q_{\mathrm{th}}`$ values for gravitational runaways due to (nonaxisymmetric) swing amplification in finite-thickness disks, as we did in Paper I for razor-thin disks; there we found $`Q_{\mathrm{th}}=`$1.3, 1.4, and 1.2 for $`\beta =\mathrm{}`$, 10, and 1 cases, respectively. If we naively assume that $`Q_c`$ is lower by $`25\%`$ in nonaxisymmetric cases as well, then thick disks would have $`Q_{\mathrm{th}}0.91`$; our results from §4 would be consistent with such thresholds. The reduction of self-gravity in finite-thickness disks also implies that stable quasi-axisymmetric spiral shock solutions can exist at a stronger spiral arm potential than predicted in the razor-thin disks (cf, Paper II).

## 4 Three-Dimensional Simulations

To investigate local three-dimensional evolution of galactic gas disks under various conditions, we have performed a total of 10 numerical simulations. Table 1 lists the important parameters for each model. Column (1) labels each run. Column (2) gives the Toomre parameter $`Q`$ (see eq. ). Note that $`Q<1`$ in a thick disk does not necessarily mean axisymmetric gravitational instability (see §3 and Fig. 3). The relative strength of self-gravity compared to external gravity is given in terms of $`s_0`$ (see eq. ) in column (3); $`s_0=25`$ models strong spiral arm regions, while $`s_0=1`$ characterizes average disk conditions. Column (4) lists the plasma parameter $`\beta `$ characterizing the magnetic field strength (see eq. ). The shear parameter $`q`$ (see eq. ) is given in column (5). Column (6) indicates whether self-gravity is implemented in the numerical evolution. Notice that even for the cases with “no” self-gravity for the evolution, the vertical total gravity profile from the “initial” equilibrium is included in order to keep the unperturbed equilibrium density profile intact. Columns (7) and (8) are our numerical resolution and the physical size of the simulation box in units of $`H`$. The relationship of $`H`$, $`Q`$, and $`s_0`$ to physical values is given in §2.2. Finally, column (9) describes the relevant physics in each model.

### 4.1 The Parker Instability Without Self-Gravity

In this section, we describe the nonlinear evolution of the Parker instability444Precisely speaking, the Parker instability that has the wavenumber vector $`𝐤`$ parallel to the initial magnetic field $`𝐁_\mathrm{𝟎}`$ is an undular mode of a more general class of magnetic buoyancy instabilities. The other branch of magnetic buoyancy instability is an interchange mode whose wavenumber vector is perpendicular to $`𝐁_\mathrm{𝟎}`$ (e.g., Hughes & Cattaneo 1987). Since the interchange mode requires $`\gamma <11/(2\beta )`$ to be unstable, where the $`z`$-independent plasma $`\beta `$ is defined by eq. (9) and $`\gamma `$ is an adiabatic index (cf, J. Kim et al. 1998), our isothermal MHD models are all stable to the interchange mode. in the absence of self-gravity. These non-self-gravitating models demonstrate the effects of rotation and shear on the Parker instability, and will later be compared with the results of self-gravitating models. From the outset, we should mention that magnetorotational instability (MRI) in shearing disks with the present parameters is not important compared to the Parker instability. It is well known that MRI of toroidal magnetic fields grows much slower than MRI of poloidal fields (cf, Balbus & Hawley 1998; Kim & Ostriker 2000). Even in a stratified disk with $`q=1.5`$ and $`\beta =12.5`$, toroidal MRI takes about ten orbital times to fully grow (Stone et al., 1996). Since the present models have weaker shear and stronger magnetic fields and are evolved only up to $`(45)t_{\mathrm{orb}}`$, effects of MRI are negligible.

#### 4.1.1 The Parker Instability in Uniformly Rotating Disks

We begin by considering the Parker instability under rigid-body rotation (model A). Figure 4$`a`$ shows the evolutionary history of the density maximum as a solid line. A three-dimensional visualization of the perturbed density with representative magnetic field lines is shown in the left panel of Figure 5. Snapshots of total surface density as well as $`x`$$`y`$ density slices are shown in Figure 6.

Initially, density perturbations which are a superposition of numerous modes with different $`𝐤`$ adjust the disk by launching MHD waves. During this relaxation process, fluid motions are very active at high-$`|z|`$ regions where the density is low. Once the system finds the most unstable Parker mode, over-dense regions slide toward the midplane along the undulating field lines, while under-dense regions move upward at the expense of gravitational potential energy (cf, Hughes & Cattaneo 1987), causing the field lines to bend further. Using the evolution of vertical velocity dispersion, we find the growth rate of the most unstable mode is $`0.4c_s/H`$ in model A. After the maximum density $`\rho _{\mathrm{max}}=1.75\rho _{00}`$ is reached at $`t/t_{\mathrm{orb}}1.8`$, the Parker instability stops growing and saturates; the corresponding maximum surface density is $`\mathrm{\Sigma }_{\mathrm{max}}=1.35\mathrm{\Sigma }_0`$. The primary reason for this is of course that the Parker instability is not a runaway process, but is eventually stabilized by tension forces from bent magnetic field lines (cf, Mouschovias 1974).

As the left frame of Figure 5 and Figure 6$`a`$ show, the dominant growing mode turns out to be antisymmetric across the midplane, and trailing with $`k_x=9(2\pi /L_x)=3.3/H`$ and $`k_y=2(2\pi /L_y)=0.74/H`$; the corresponding wavelengths are $`\lambda _x=2H`$ and $`\lambda _y=8.5H`$. The antisymmetric (in density) property of the dominant Parker mode in our model disks is consistent with the results of Santillán et al. (2000), who showed that numerical simulations with random initial perturbations generally exhibit antisymmetric modes, even though the linear growth rates for symmetric and antisymmetric modes of the Parker instability are almost indistinguishable in galactic environments (J. Kim et al., 2000).555Horiuchi et al. (1988) showed that the growth rates of antisymmetric and symmetric modes depend on their $`ϵ`$ parameter (see eq.  of Horiuchi et al. 1988). In galactic conditions, however, $`ϵ10^3(R_0/10\mathrm{kpc})^2(\mathrm{\Omega }/26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1)^2(c_s/7\mathrm{km}\mathrm{s}^1)^2`$, so that their Fig. 5 reads that the maximum growth rates of antisymmetric modes are almost the same as those of symmetric modes.

The dominant mode in model A is trailing ($`k_x/k_y>0`$), as well. Linear dispersion relations for the Parker instability with rotation are independent of the sign of $`k_x/k_y`$ (cf, Shu 1974), and we have found that for low-resolution simulations seeded with either leading or trailing modes, growth rates are similar. However, for the random perturbation form we adopt, trailing modes grow sooner because of a shorter readjustment phase, for the rotating cases. Additional tests for the nonrotating case show that the growth of trailing and leading modes is indistinguishable.

Three-dimensional Parker instability attains its maximum growth rate at $`k_x\mathrm{}`$, because larger $`k_x`$ provides more space for a raised section of a flux tube to expand in the radial direction (Parker, 1967). As pointed out by Shu (1974), however, the growth rate of $`k_x=\mathrm{}`$ mode is not very different from that of the $`k_x=0`$ modes and the growth rate is essentially constant for $`k_xH\stackrel{>}{}5`$. Since our high resolution simulations can resolve modes up to $`k_xH12`$, the dominance of the $`k_xH=3.3`$ modes in model A is physical rather than numerical. The dominance of this particular mode over potential rivals may in part be attributed to the shape of initial perturbation spectrum which has larger power at smaller $`k_x`$, while larger-$`k_x`$ modes, if resolved, grow faster; a compromise may then select intermediate-$`k_x`$ modes that have large enough initial power and fast enough growth rates. It is also possible that nonlinear effects help to select this mode, perhaps preferring $`\lambda _x`$ close to the full disk thickness. We note that although uniform rotation tends to stabilize the Parker instability since the Coriolis force couples azimuthal and radial motions, its effect is not significant if $`k_x/k_y1`$ (cf, Shu 1974; Zweibel & Kulsrud 1975; Foglizzo & Tagger 1994).

In addition to helping select trailing modes for initial growth, the effects of rotation are evident in other details of the model evolution. As material falls toward the midplane along wavy magnetic field lines, it also gains speed in the $`\widehat{𝐲}`$ direction, and its path is consequently deflected in the $`\widehat{𝐱}`$ direction by the Coriolis force. Gas moving toward a magnetic valley with positive (negative) $`v_y`$ is redirected radially outward (inward). The net result is that over-dense regions at valleys tend to rotate in the counterclockwise direction, which in turn twists magnetic field lines in the valleys (cf, Kim, Ryu, & Jones 2001). At the point when the instability begins to saturate, gas at the valleys is slightly overcompressed and begins to reexpand vertically. At the same time, the fields lines that previously became twisted due to the Coriolis forces now exert back torques, inducing gas motions in the clockwise direction. The meridional flows at the disk midplane, as displayed in Figure 6$`c`$, show this sort of behavior. The flows are generally subsonic near the midplane, while high altitude regions ($`|z|/H>2.5`$) have supersonic velocities.

Evolution of the system following saturation is rather complicated, involving nonlinear interactions of numerous MHD waves. Overpressured, dense regions shown in Figure 6$`b`$ soon expand in the horizontal direction. Expanding material interacts with gas that moves along the vertical direction. These active, but still mostly subsonic, flow motions eventually homogenize the density, resulting in a smooth surface density distribution (Fig. 6$`d`$). The maximum surface density after 3 orbits is only 1.13 times the initial value.

The small-scale chaotic flow associated with nonlinear development of large-$`|𝐤|`$ modes is also likely to cause reconnection of the curved and twisted magnetic field lines, preferentially at low-$`|z|`$ regions (cf, J. Kim et al. 1998; Kim, Ryu, & Jones 2001). This in turn reduces the field strength near the midplane, which was already decreased during the process of (midplane-crossing) Parker instability.

We define the differential mass

$$\frac{dM(z)}{dz}\rho 𝑑x𝑑y,$$
(19)
and the differential magnetic flux averaged along the $`y`$-direction

$$\frac{d\mathrm{\Phi }_B(z)}{dz}\frac{1}{L_y}B_y𝑑x𝑑y,$$
(20)
and plot $`dM/dz`$, $`d\mathrm{\Phi }_B/dz`$, and $`dM/d\mathrm{\Phi }_B`$ as functions of $`z`$ in Figure 7 for five time epochs. For $`t/t_{\mathrm{orb}}\stackrel{<}{}1`$, both profiles remain unchanged. As the Parker instability develops considerably ($`1\stackrel{<}{}t/t_{\mathrm{orb}}\stackrel{<}{}2`$), field lines float upwards, supplying magnetic flux at high altitude regions; magnetic fields may escape the simulation box because of the outflow $`z`$-boundary conditions we adopted. Magnetic fields rise up and gas slips down, enhancing the mass-to-flux ratio further at the midplane (Fig. 7$`c`$), and evening out the $`z`$-distribution of magnetic flux (Fig. 7$`b`$), but mass enhancement near the midplane due to the Parker instability is generally quite moderate (Fig. 7$`a`$). At the end of the simulation, the overall density profile in model A is very close to another hydrostatic equilibrium where the initial effective external force ($`g_{\mathrm{eff}}`$; see eq. ) counterbalances the thermal pressure gradient alone. As Figure 7$`b`$ suggests, the mean magnetic field (which is predominantly azimuthal) is nearly constant with height, and therefore does not contribute much in the final force balance. At $`t/t_{\mathrm{orb}}=4`$, the density-weighted velocity dispersions are calculated to be $`\sigma _x\rho v_x^2^{1/2}/\rho ^{1/2}2.5\times 10^2c_s`$, $`\sigma _y\rho v_y^2^{1/2}/\rho ^{1/2}7.3\times 10^2c_s`$, and $`\sigma _z\rho v_z^2^{1/2}/\rho ^{1/2}2.4\times 10^2c_s`$, where $``$ denotes a volume average. Given the low level of the velocity dispersion, the Parker instability is unlikely a significant source of interstellar turbulence in galactic disks.

#### 4.1.2 Effects of Differential Rotation

One of the crucial impacts differential rotation makes on the evolution of the Parker instability is that it changes the power spectrum of initial perturbations. In a shearing disk, the wavefront of any disturbance evolves in response to the kinematics of the background flow; for uniform shear with $`qd\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}R`$, the radial wavenumber of any pattern linearly increases with time as $`k_x=k_x(0)+q\mathrm{\Omega }k_yt`$. In our chosen initial density perturbation spectrum, the mode with $`k_x(0)=0`$ and $`k_y=1(2\pi /L_x)=0.37/H`$ has much larger power than higher-$`k_x`$ modes at the same $`k_y`$. The higher-$`k_x`$ modes, however, have higher growth rates than the $`k_x=0`$ modes. As a consequence of the initial power distribution, non-shearing model A thus has to spend about one orbital time before the (large-$`k_x`$) Parker modes grow enough to produce significant density fluctuation (Figure 4$`a`$). On the other hand, the background shear in model B can transform the initial $`k_x(0)=0`$ mode, without changing its power, into a higher-$`k_x`$ disturbance – for which the Parker instability is more efficient. Therefore, in Figure 4$`a`$ the initial growth of perturbations is more evident in model B than in model A. Structural analysis shows that this initially growing mode in model B is a $`z`$-symmetric (in density) mode with $`k_yH=0.37`$ ($`\lambda _y=17H`$) and confined to $`|z|/H<1`$. In contrast, the most unstable mode in model A is a $`z`$-antisymmetric mode with $`k_yH=0.74`$. The growth rate of the $`k_yH=0.37`$ mode in model B is found to be $`0.2c_s/H`$, which is just half of the growth rate of model A.

While the initial evolution of shearing model B is governed by the growth of $`z`$-symmetric modes near the midplane, $`z`$-antisymmetric perturbations that have lower initial power (since $`k_z`$ is larger) also grow in high altitude regions at double the rate. As they grow, over-dense regions fall down toward the midplane and begin to interact with the symmetric modes. The right panel of Figure 5 displays perturbed density and magnetic field lines in model B at $`t=10\mathrm{\Omega }^1`$. The condensations evident at the midplane are associated with the symmetric mode with $`k_x=10(2\pi /L_x)`$ ($`\lambda _x=2H`$). Below or above these are “tadpole-shaped” condensations associated with the antisymmetric mode having $`k_x=20(2\pi /L_x)`$ ($`\lambda _x=H`$). Magnetic field lines from model B appear straighter than from model A, suggesting that the (longer-wavelength) symmetric modes partially prevent the (shorter-wavelength) antisymmetric modes from crossing the midplane. As antisymmetric modes grow dominant later on ($`t/t_{\mathrm{orb}}>1.7`$), field lines become progressively more curved, but never to the extent seen in model A.

For shearing model B, antisymmetric modes culminate their growth at $`t/t_{\mathrm{orb}}2.1`$. After model saturation, the system experiences vigorous interactions between symmetric and antisymmetric modes. An added complication is the background shear, which stretches out any condensation that forms. We plot snapshots of surface density $`\mathrm{\Sigma }/\mathrm{\Sigma }_0`$ in Figure 8, which shows that at $`t/t_{\mathrm{orb}}=2`$ the system is dominated by high-$`k_x`$ modes, while surface density at $`t/t_{\mathrm{orb}}=3`$ is rather uniformly distributed. To some extent, the smooth surface density is attributable to our limited numerical resolution. However, even if the true spatial scale at which dissipative processes smooth out strongly phase-mixed structures is small compared to our numerical resolution, there are evidently no nonlinear feedback processes to drive late-time growth of large-scale perturbations. Overall, the late time character of the system is a new, large-scale vertical equilibrium containing small-scale, fluctuating density and velocity fields. The density-weighted dispersions of the perturbed velocity in model B at time $`t/t_{\mathrm{orb}}=4`$ are $`\sigma _x0.10c_s`$, $`\sigma _y0.15c_s`$, and $`\sigma _z0.07c_s`$, which are too small to be of much significance for the observable interstellar turbulence in galactic disks.

As mentioned before, symmetric modes in the shearing model B impede the midplane penetration of antisymmetric modes, resulting in less inflated field lines. The antisymmetric modes have half the azimuthal wavelength of the symmetric modes, and also shorter radial wavelengths, which helps perturbations squeeze across the midplane. Nevertheless, the escape of magnetic flux into high altitude regions is significantly reduced. Figure 9 shows the variations of $`dM/dz`$, $`d\mathrm{\Phi }_B/dz`$, and $`dM/d\mathrm{\Phi }_B`$ for model B. Compared with Figure 7 for non-sheared model A, Figure 9 for model B illustrates that shear in the Parker instability results in a configuration where magnetic flux is enhanced in two layers ($`0.5H\stackrel{<}{}|z|\stackrel{<}{}1.5H`$) surrounding the midplane, while the corona has lower $`B`$ (higher $`\beta `$) than the non-sheared case (cf, Shibata, Tajima, & Matsumoto 1990).666On the other hand, Miller & Stone (2000) showed that magnetorotational instability in a shearing disk tends to produce a strongly magnetized corona. Apart from the shearing background velocity, small-scale motions are chaotic, which can facilitate the reconnection of twisted magnetic field lines. The “layered” magnetic profile may in part arise from enhanced reconnection near the midplane, raising $`dM/d\mathrm{\Phi }_B`$ at $`|z|/H\stackrel{<}{}0.5`$ (Fig. 9$`c`$).

To summarize, we do not find any firm evidence that shear significantly stabilizes the Parker instability. This may be because Parker modes are the most unstable at large $`k_x`$. By kinematically shifting the $`x`$-wavenumbers of perturbations at lowest $`k_z`$ and $`k_y`$ which have the largest initial amplitudes – and are associated with symmetric modes – the introduction of shear gives an early boost to symmetric modes. The more slowly-growing symmetric modes are later dominated by shorter-wavelength antisymmetric modes that grow more strongly (but start from lower amplitudes). The midplane-centered symmetric modes nevertheless interfere with the midplane-crossing antisymmetric modes, reducing the latter’s amplification and limiting the inflation of field lines toward high latitudes. Except for having slightly more (less) magnetic support at moderate (low) $`|z|`$, the final state of the shearing model does not differ much overall from that of the nonshearing model, consisting of a new vertical equilibrium with enhanced central concentration.

### 4.2 Gravitational Instability in Shearing Thick Disks

#### 4.2.1 Swing Amplification in Unmagnetized Disk

Self-gravitating evolution of an unmagnetized, thick disk is generally very similar to that of a razor-thin disk, because without magnetic fields the induced motions are almost planar. One qualitative difference the disk’s finite thickness makes is that it tends to be stabilizing by reducing self-gravity at the midplane (see §3).

As shown in Figure 4$`b`$ as a dotted line, the initial evolution of density in model C is governed by nonaxisymmetric swing amplification. With $`Q=0.7`$ and $`s_0=25`$, hydrodynamic model C is stable to axisymmetric perturbations (Fig. 3). As wavefronts swing around from leading to trailing, epicyclic motion with the same rotational sense allows overdense fluid elements to linger in wave crests. This extended exposure to self-gravity enables perturbations to grow until the radial wavenumber becomes large (Goldreich & Lynden-Bell 1965b; Toomre 1981; see also Paper I). The swing amplification in model C saturates at $`t/t_{\mathrm{orb}}0.8`$, producing a state of nonlinearly interacting sheared wavelets of relatively high $`k_x`$. Although the medium continues to shear kinematically, nonlinear interactions maintain a dominant population of wavelets at moderate $`k_x`$. As displayed in the left panel of Figure 10, for example, the surface density map at $`t/t_{\mathrm{orb}}=2`$ has the most power at $`k_x=7(2\pi /L_x)`$.777Fourier analysis shows that modes of $`k_x>10(2\pi /L_x`$) have negligible power throughout the saturation stage, confirming that our simulation results are not limited by numerical resolution. At this time, $`|\delta v_y|0.5|v_x|`$ and $`|v_z|0.04|v_x|`$, indicating the vertical motions are indeed insignificant compared with the planar motions.

Much as we found for two-dimensional models (Paper I), the evolution of model C after saturation indicates that nonlinear interactions repopulate small-$`|k_x|`$ modes, which then rejuvenate swing amplification. This second growth phase yields a few high-density filaments (middle panel of Fig. 10). The filaments collide with each other and become gravitationally unstable, eventually forming two gravitationally bound clumps at the end of simulation. The mass of each clump is $`7\%`$ of the total mass, corresponding to $`4\times 10^6\mathrm{M}_{}`$ (right panel of Fig. 10); this is twice the corresponding two-dimensional Jeans mass, $`M_\mathrm{J}c_s^4/(G^2\mathrm{\Sigma }_0)=2\times 10^6\mathrm{M}_{}(c_s/7.0\mathrm{km}\mathrm{s}^1)^4(\mathrm{\Sigma }_0/65\mathrm{M}_{}\mathrm{pc}^2)^1`$, consistent with the results of Paper I where we found the fragments formed in a two-dimensional disk model with comparable $`Q`$ and surface density typically have $`(1/22)M_\mathrm{J}`$ (Model H06 in Paper I).

In Paper I, we reported that when $`Q`$ is sufficiently small, razor-thin disks form gravitationally bound clumps via three different types of secondary instabilities. These include parallel fragmentation ($`Q\stackrel{<}{}0.8`$), gravitationally induced collisions of nonlinear sheared wavelets ($`0.8\stackrel{<}{}Q\stackrel{<}{}1.1`$), and rejuvenated swing amplification ($`Q\stackrel{>}{}1.1`$ for unmagnetized disks). Had model C been evolved in the razor-thin approximation, it would have experienced gravitational runaway via parallel fragmentation long before initial swing amplification saturated. It is remarkable that a thick disk with $`Q`$ as small as 0.7 is barely unstable, requiring two successive phases of swing amplification before runaway occurs.

#### 4.2.2 Parker-Swing Instability

The models of §4.1.2 and §4.2.1 represent systems where either gravity or magnetic fields are entirely turned off, to investigate the physics of sheared Parker and thick-disk swing instabilities separately. Our primary interest here, however, is to understand how magnetic buoyancy may interact with the self-gravitating instabilities already known to exist from two-dimensional studies (Paper I). Model D focuses on the swing-Parker interaction, incorporating self-gravity, shear, and thermal-equipartition toroidal magnetic fields. Time evolution of the maximum density in model D is plotted in Figure 4$`b`$ as a solid line. Clearly, the initial phase of evolution is dominated by swing amplification, with very similar behavior to model C (dotted line). Unlike in models A and B of the Parker instability without self-gravity, density in model D begins to grow immediately as $`k_x`$ increases. Embedded magnetic fields reduce the amplification factor in model D by 63% compared to model C.

Swing amplification saturates at about $`t/t_{\mathrm{orb}}0.8`$ when $`k_x`$ becomes sufficiently large. With embedded azimuthal magnetic fields, however, model D is also subject to the Parker instability. The shearing wavelets grown from the initial swing amplification in fact strengthen the symmetric Parker modes. Density in model D increases steadily for $`1\stackrel{<}{}t/t_{\mathrm{orb}}\stackrel{<}{}3`$; model B shows similar behavior at earlier times in Figure 4$`a`$. At the same time nonlinear interaction of the wavelets – as discussed in §4.2.1 – keeps the power spectrum centered at moderate $`k_x`$. Figure 11 shows surface density snapshots at $`t/t_{\mathrm{orb}}=2`$, 3, and 4.1 of model D. At $`t/t_{\mathrm{orb}}=2`$, most of the power is concentrated in the mode with $`k_x=9(2\pi /L_x)`$ and $`k_y=1(2\pi /L_x)`$. Vertical flows due to the Parker modes couple with in-plane (primarily radial) motions, enabling more vigorous nonlinear feedback. Rejuvenated swing amplification of small-$`|k_x|`$ modes follows, and a few dense filaments form (middle of Fig. 11). They later experience mutual collisions, ultimately forming two gravitationally bound condensations. (right of Fig. 11). Each clump has 8% (left) and 12% (right) of the total mass, corresponding approximately to $`9.6\times 10^6\mathrm{M}_{}`$ and $`1.4\times 10^7\mathrm{M}_{}`$.

Because of the strong $`z`$-symmetric perturbations resulting from the initial swing amplification, model D does not provide fertile ground for growth of antisymmetric Parker modes. The latter may grow slightly in high- and middle-altitude regions, but never become dynamically important. Consequently, magnetic field lines do not float upward as much as in the non-self-gravitating model B. We draw $`dM/dz`$, $`d\mathrm{\Phi }_B/dz`$, and $`dM/d\mathrm{\Phi }_B`$ for model D in Figure 12, which illustrates that differential mass and flux at $`|z|/H\stackrel{>}{}1`$ are almost unchanged for $`t/t_{\mathrm{orb}}\stackrel{<}{}2`$. Nonlinear interaction of wavelets may cause magnetic field lines to reconnect, reducing $`d\mathrm{\Phi }_B/dz`$ near the midplane.

Model D investigated how Parker instability interacts with swing amplification under conditions of relatively high gaseous surface density such that the (height-integrated) $`Q`$ parameter is 0.7 and gas gravity dominates stellar gravity locally ($`s_0=25`$). Regions in which both gas surface density and shear are so high may, however, be exceptional in contemporary galaxies. To investigate the effects of the Parker-swing interaction under other galactic conditions, we also perform low-resolution simulations for two cases with lower gaseous surface density. Both of these models have $`s_0=1`$; model G has $`Q=1.5`$ and model I has $`Q=1.0`$. For each of these, we also perform simulations in which self-gravity is turned off (models F and H, respectively) to provide “pure Parker” comparison models. In addition, we turn off magnetic fields in model J to provide a $`Q=1.0`$ “pure swing” counterpart to model I. The resulting evolutionary histories of maximum density in models F–J are presented in Figure 13.

When $`Q=1.5`$ and $`s_0=1`$, corresponding approximately to average disk conditions, self-gravity is so weak that the growth of density via swing amplification in model G (solid line in Fig. 13$`a`$) is by less than a factor of two. Since the amplitude of the perturbed density field remains small, nonlinear feedback from wave interactions do not inhibit the kinematic increase of $`k_x`$. The medium keeps shearing out, and Parker modes slowly emerge during the later part of evolution ($`t/t_{\mathrm{orb}}\stackrel{>}{}2.5`$). Model F, in which self-gravity is turned off (dotted line in Figure 13$`a`$), has late-stage evolution after saturation of Parker modes very similar to that of model G. Namely, a new quasi-static vertical equilibrium, with low-amplitude velocity and density fluctuations, is achieved. This state is similar to that of late-stage model B, except that the fluctuation level is lower by a factor of three. The reason for the differences between the two non-self-gravitating models (B and F) is that the former has a smaller scale height due to a stronger gas self-gravity in the initial equilibrium state888As noted in the beginning of §4, our non-self-gravitating (for evolution) models still have the initial equilibrium self-gravity that together with external gravity balances thermal and magnetic pressure forces in the background state. Although the treatment of self-gravity is not self-consistent in model B, its primary function is as a control for the otherwise identical model D, in order to assess the small-scale evolutionary effects produced by self-gravity compared to Parker modes alone. ($`H\mathrm{\Omega }/c_s`$ curves in Fig. 1), thereby enabling the Parker instability to grow faster and more strongly (cf. eq. ). This is also why the clear transition from symmetric to antisymmetric Parker modes that occurs at $`t/t_{\mathrm{orb}}1.7`$ in model B (Fig. 4$`a`$, dotted line) is not obvious in the density history of model F (Fig. 13$`a`$, dotted line). The conclusive result from the $`Q=1.5`$, $`s_0=1`$ model is that neither the Parker instability nor self-gravity could lead to significant density condensation under “average disk” conditions that are dominated by equipartition-amplitude, large-scale toroidal magnetic fields.

Models H, I, and J investigate the Parker/swing connection under conditions of higher disk surface density such that $`Q=1.0`$; we set $`s_0=1`$. The evolutionary histories of $`\rho _{\mathrm{max}}/\rho _{00}`$ are shown in Figure 13$`b`$. Evidently, comparison of the development of magnetized (solid line; model I) and unmagnetized (dashed line; model J) models shows that Parker modes can help push a region into self-gravitational runaway when $`Q`$ is at a “marginal” level. Model I produces a couple of dense filaments at $`t/t_{\mathrm{orb}}=4.8`$ and is expected to end in collapse, whereas the initial swing amplification in model J phase-mixes away within $`3`$ orbits. It is evidently a combination of the extra self-gravity produced in “valleys” of Parker modes combined with the vorticity-removing ability of magnetic fields that tips the balance in allowing a rejuvenated, runaway swing phase in model I. Model H, in which self-gravity is turned off (dotted line), shows the development of density perturbations for the “pure Parker” case, which is quite similar to that of model B (see Fig. 4$`a`$), even though $`s_0`$ and $`Q`$ differ.

In conclusion, the models of this section show that the Parker instability does not significantly alter the ability of a large-scale region of a galactic disk to form gravitationally bound clouds. Under “average disk” conditions ($`Q=1.5`$, $`s_0=1`$, $`q=1`$), the disk turns out to be extremely stable, with only mildly fluctuating density in the final state. Both magnetized and unmagnetized disks can form bound clouds for the high gas surface-density conditions ($`Q=0.7`$, $`s_0=25`$) that represent spiral arm regions (except using “normal” shear). For intermediate surface density ($`Q=1.0`$), only magnetized disks form bound clouds. While we have not pinpointed a threshold $`Q`$ value for combined Parker-swing nonlinear runaway, models G (see Fig. 13$`a`$) and I (see Fig. 13$`b`$) taken together bracket its value between 1.5 and 1.0 (when $`s_0=1`$). Even when bound clouds are created, the Parker instability plays only an auxiliary role by helping secondary swing amplification; the formation time scale and spacing of clouds are quite different from what the linear-theory Parker analysis predicts.

### 4.3 Magneto-Jeans Instability

The final issue we address in this paper is the development of the magneto-Jeans instability (MJI) in a three-dimensional disk. We previously studied MJI under uniform weak-shear conditions (Paper I), and under self-consistent spiral arm conditions (Paper II), both in a height-integrated two-dimensional idealization. While the latter varying shear and varying density approach is more realistic, in this first three-dimensional study we adopt the simplified uniform model, which is adequate to capture the basic physical processes: model E has the same high-gas-density-parameter set ($`Q=0.7`$, $`s_0=25`$) used in model D, but sets $`q=0`$ to represent better the reduced (or negative) shear that prevails within spiral arms (see Paper II).

Physically, MJI is a form of nonaxisymmetric Jeans instability in rotating disks containing embedded azimuthal magnetic fields. Transverse (Alfvénic) magnetic perturbations produce radial tension forces that resist the stabilizing effect of the Coriolis force; magnetic transfer of angular momentum allows secular contraction instead of epicyclic motion (Lynden-Bell 1966; Elmegreen 1987; Paper I). Although MJI is similar to the Parker instability in that both rely on background azimuthal magnetic fields and nonaxisymmetric perturbations, they are otherwise completely different. The condition of low shear is essential for MJI, because magnetosonic waves begin to quench the Jeans instability when kinematic shear increases $`k_x`$; the Parker instability, in contrast, prefers high $`k_x`$. Vertical motions are key in the Parker instability, whereas MJI involves primarily planar motions. Since both have been proposed as possible ways to trigger formation of giant clouds in spiral arms, it is interesting to investigate nonlinear evolution when both instabilities can occur.

Based on model E, we conclude that MJI, if present, always dominates the Parker instability. Although equation (12) suggests the Parker instability can have a high initial growth rate, its development never produces highly nonlinear near-midplane density perturbations (see solid line in Fig. 4$`a`$). As the dashed line in Figure 4$`b`$ shows, the density in model E ($`Q=0.7`$) starts to grow immediately and this growth never saturates, just as we found for MJI in Paper I. We have also performed low-resolution simulations for MJI in nonshearing disks with $`Q=1.0`$ and 1.5 (not listed in Table 1), and confirmed that MJI can drive them into gravitational runway. The time scale for bound cloud formation becomes longer as $`Q`$ increases, so that disks with $`Q=1.0`$ and 1.5 form dense clumps within 1 and 2 orbits, respectively, while model E with $`Q=0.7`$ takes less than one orbital time to become fully unstable (Fig. 4$`b`$).

In Figure 14, we plot perspective views of isodensity surfaces and magnetic field lines, together with midplane distributions of density, velocity, and magnetic fields, for model E at $`t/t_{\mathrm{orb}}=0.5`$. As the perspective view shows, the magnetic fields run almost parallel to the $`x`$$`y`$ plane near $`z=0`$, as do the velocity vectors even during highly nonlinear stages. Centered at $`(x/L_x,y/L_y,z/L_z)=(0.10,0.19,0)`$ and $`(0.14,0.05,0)`$, two condensations are produced within one orbital time; these each contain 27% and 20% of the total mass, corresponding roughly to $`3\times 10^7\mathrm{M}_{}`$.

As discussed in §3.2 for axisymmetric gravitational instabilities, the dilution of in-plane self-gravity due to the finite scale height also affects the most unstable scale lengths and growth rates of MJI. The simple replacement $`\mathrm{\Sigma }_0k2\rho _{00}Hk/(1+Hk)`$ is sufficient to take into account this geometrical dilution effect in an approximate fashion. Utilizing equation (23) of Paper I (see also Lynden-Bell 1966), the instantaneous dispersion relation for MJI in thick disks becomes

$$0=\omega ^4\left[\kappa ^2+(c_s^2+v_\mathrm{A}^2)k^24\pi G\rho _{00}\frac{kH}{1+kH}\right]\omega ^2+v_\mathrm{A}^2k_y^2\left(c_s^2k^24\pi G\rho _{00}\frac{kH}{1+kH}\right),$$
(21)
where $`k^2=k_x(t)^2+k_y^2`$ and $`k_x(t)=k_x(0)+q\mathrm{\Omega }k_yt`$. Perturbations require $`k_y0`$ (nonaxisymmetric) and $`kH_{\mathrm{}}(1+kH)<2`$, or $`k<k_c(\sqrt{1+8H/H_{\mathrm{}}}1)/2H`$ to be unstable to MJI (see eq.  for the definition of $`H_{\mathrm{}}`$). Note that the critical wavenumber $`k_c`$ of MJI is independent of magnetic pressure (except through $`H`$) and rotation, so that MJI operates at smaller scales along the magnetic field than axisymmetric gravitational instability. Since $`H/H_{\mathrm{}}1`$ for $`s_0\stackrel{>}{}1`$ (Fig. 1), the condition for MJI is roughly $`kH\stackrel{<}{}1`$, which is similar to the instability criterion for the Parker modes.999In investigating the formation of filamentary structures in the ISM, Chou et al. (2000) studied disk models that are similar to model E, that is, magnetized, uniformly-rotating, and self-gravitating. They identified the structure forming mechanism with the Parker-Jeans instability, but we believe that while the Parker mechanism may contribute during the linear growth phase, the primary driver of nonlinear structure formation is the MJI.

Using equation (21) and the results of §3.1, we obtain the peak ($`k_x=0`$) growth rates for MJI for model E ($`Q=0.7`$, $`s_0=25`$ and $`\beta =1`$) as a function of $`k_y`$, and plot them in Figure 15. For comparison, we also plot the growth rates for the corresponding razor-thin disk. The maximum growth rate for a thick disk is $`1.12\mathrm{\Omega }`$ at $`k_{y,\mathrm{max}}=0.74/H`$, so $`\lambda _{y,\mathrm{max}}8H`$. Since $`L_y17H`$ for our box, this would predict $`\lambda _{y,\mathrm{max}}/L_y0.5`$, i.e., two $`x`$-aligned filaments. This is indeed consistent with the density structure in Fig. 14, although there is also a small contribution from the $`\lambda _{y,\mathrm{max}}=L_y`$ mode that has higher power from our initial conditions. Note that the most unstable $`\lambda _y`$ of MJI is in fact the same as that of the dominant Parker mode that grows in non-sheared model A, although that mode is midplane-antisymmetric in density. For both $`x`$\- and $`y`$-scales comparable to $`\lambda _{y,\mathrm{max}}`$, the corresponding mass is $`M=\mathrm{\Sigma }_0(2\pi /k_{y,\mathrm{max}})^23\times 10^7\mathrm{M}_{}(\mathrm{\Sigma }_0/65\mathrm{M}_{}\mathrm{pc}^2)(H/80\mathrm{pc})^2`$, consistent with the masses of clumps formed in model E.101010Because $`H`$ varies inversely with $`\mathrm{\Sigma }`$ from eq. (15), the condensation mass would be proportional to $`\mathrm{\Sigma }^1`$. From Figure 15, we note that a razor-thin disk would characteristically respond to MJI twice as fast, at half the length scale, of a thick disk. Consequently, the characteristic Jeans mass is about four times bigger in a thick disk than the thin-disk prediction. We note that this result is consistent with simple prediction for an unmagnetized disk using the vertical equilibrium condition of equation (15): the central density is $`\pi G\mathrm{\Sigma }_0^2/2c_s^2`$, and with the three-dimensional Jeans wavelength $`\lambda _{\mathrm{J},3\mathrm{D}}=c_s(\pi /G\rho )^{1/2}`$, the nominal three-dimensional Jeans mass would be $`M_{\mathrm{J},3\mathrm{D}}=\pi \sqrt{2}c_s^4/(G^2\mathrm{\Sigma }_0)`$. With $`\lambda _{\mathrm{J},2\mathrm{D}}=c_s^2/G\mathrm{\Sigma }_0`$, the two-dimensional Jeans mass is $`M_{\mathrm{J},2\mathrm{D}}=c_s^4/(G^2\mathrm{\Sigma }_0)0.23M_{\mathrm{J},3\mathrm{D}}`$.

## 5 Conclusions

### 5.1 Summary

A variety of evidence supports the idea that giant, star-forming clouds form as the result of large-scale collective effects – i.e., dynamical instabilities. Candidate mechanisms that have been proposed on the basis of linear-theory analysis include the swing amplifier (Goldreich & Lynden-Bell, 1965b; Julian & Toomre, 1966), magneto-Jeans instabilities (Lynden-Bell 1966; Elmegreen 1987; Paper I), and the Parker instability (Parker, 1966). The first two of these involve primarily in-plane motions and are fundamentally driven by self-gravity, with either shear or magnetic tension acting to neutralize the stabilizing tendency of the Coriolis force. For the third mechanism, vertical motion, driven by magnetic buoyancy, is crucial.

Based on linear-theory growth rates and characteristic spatial scales, all of the above mechanisms could in principle be dynamically important under various galactic conditions. In practice, because giant clouds represent highly overdense regions compared to mean ISM properties, it is important to understand how the candidate instability mechanisms develop in the nonlinear regime: i.e., whether and/or how an instability saturates. Addressing these questions requires direct numerical simulations to evolve the fluid dynamics equations into the nonlinear regime.

In previous work, we have investigated the nonlinear development of the swing and magneto-Jeans mechanisms in a thin-disk (two-dimensional) approximation (Papers I, II). Here, we extend those studies by allowing for fully three-dimensional dynamics, which also enables us to investigate the nonlinear development of the Parker instability allowing for rotation, shear, and self-gravity. Our primary goals were to compare dynamical development, and especially “final-state” outcomes, of model disks under various conditions. By controlling for different physical effects (i.e., turning gravity, rotation, shear, and magnetic fields “off” and “on”), and exploring parameter space, we were able to assess the potential consequences for bound cloud formation of the various instabilities alone and in combination.

The numerical models we adopt are radially localized, vertically stratified, shearing, isothermal disks with embedded magnetic field lines (see §2.1). The initial rate of shear in the background flow is measured by $`qd\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}R`$, and this overall shear gradient is maintained by enforcing sheared-periodic radial boundary conditions; we study cases with $`q=1`$ and 0. The initial magnetic field is assumed to be azimuthal and its vertical stratification is determined such that the Alfvén speed $`v_\mathrm{A}`$ is spatially constant. The density and magnetic field strength in our model disks are characterized by the Toomre parameter $`Q`$ (see eq. ) and the plasma parameter $`\beta `$ (see eq. ), respectively. Magnetic fields are either initiated at an equipartition level ($`v_\mathrm{A}/c_s=1`$)111111Thus magnetorotational instability cannot grow; see below., or set to zero. We set up initial magnetohydrostatic equilibria with thermal and magnetic pressure gradients balancing gaseous self-gravity and external gravity from a stellar component. We adopt a linear model for the external (stellar-disk) gravity (see eq. ) and parameterize its magnitude using $`s_0`$ (see eq. ), which corresponds roughly to the ratio of the gas disk’s vertical gravity to the external gravity at one scale height. The fiducial sets of parameters our model disks adopt are $`Q=1.5`$, $`s_0=1`$ for the “average disk”, and $`Q=0.7`$, $`s_0=25`$ for highly compressed spiral arm regions; we also perform additional simulations with $`Q=1`$ and $`s_0=1`$.

Our simulation boxes are 8 scale-heights ($`H`$) thick, with area either $`17H\times 17H`$ or $`25H\times 25H`$ in the $`x`$$`y`$ plane. The relation between dimensionless simulation parameters and physical scales is explained in §2.2. For all our simulations, we introduce a spectrum of low-amplitude random perturbations, and numerically integrate the ideal MHD equations up to four or five orbital periods. In §3, we first constructed the density profiles for vertical static equilibria and calculated the variation of disk thickness as a function of $`s_0`$ (§3.1), and then studied axisymmetric gravitational instability of these thick disk models (§3.2). Finite disk thickness reduces the critical value of the Toomre stability parameter relative to the zero-thickness limit, and also increases the spatial scale required for instability. We found that for $`s_0=1`$, $`Q_c0.75`$, 0.72, and 0.57 for $`\beta =\mathrm{}`$, 10, and 1 cases, respectively, suggesting that the average disk is highly stable to axisymmetric perturbations. We also showed that a simple modification of the axisymmetric dispersion relation to account for finite thickness (eq. ) yields excellent agreement with the results of simulations for the value of $`Q_c`$.

In §4, we present nonlinear evolution of fully three-dimensional disk models, first considering the Parker instability in isolation (§4.1), then investigating evolution of the swing and Parker-swing mechanisms (§4.2), and finally considering how self-gravitating perturbations grow under low shear conditions (§4.3).

The main results drawn from these three-dimensional numerical simulations can be summarized as follows:

1. In the absence of self-gravity, the dominant growing Parker mode under uniform rotation is found to be antisymmetric in density with respect to $`z=0`$ and tailing with radial wavelength $`\lambda _x=1.9H`$ and azimuthal wavelength $`\lambda _y=8.5H`$, and has a growth rate $`0.4c_s/H`$. Density fluctuations grow as over-dense regions slide toward the midplane along undulating magnetic field lines and unloaded portions of the field buoyantly rise, but growth slows as nonlinear saturation sets in. The maximum density at saturation is less than twice the initial midplane density. After saturation of the dominant Parker mode, nonlinearly developed small-scale modes interact with each other, dispersing density structures characteristic of the Parker instability. Escape of magnetic flux across the vertical boundaries and likely reconnection of the twisted magnetic fields prompted by small-scale motions raise the mass-to-flux ratio near the midplane. The system progressively moves to a new global vertical equilibrium in which the magnetic field is nearly uniform, hence magnetic pressure support against gravity is negligible. In the saturated state, small scale fluctuating velocities are still present, but their density-weighted dispersions are less than 10% of the isothermal speed speed.

2. When a disk is subject to differential rotation, midplane symmetric (in density) Parker modes with $`\lambda _y=17H`$ are the first to dominate the evolution. This initial dominance is explained by the large initial amplitude of longer-wavelength perturbations, combined with the ability of shear to increase $`k_x`$ to the regime where Parker instabilities are more efficient. Smaller-$`\lambda _y`$ antisymmetric modes start with lower amplitudes than the symmetric modes but grow more rapidly, moving gradually downwards from high altitudes to interact with the symmetric mode. The mode interaction reduces magnetic field line inflation and loss compared to the non-shearing case. The final state of the non-self-gravitating shearing disk model is similar to that of the non-shearing disk, consisting of an overall vertical equilibrium with small-scale, low-amplitude fluctuations in density and velocity fields.

3. When we include both self-gravity and uniform shear in our models, the nonlinear evolution of structure is principally dominated by the swing amplification mechanism, although Parker modes play a role in seeding structure. Low-$`Q`$ ($`Q=0.7`$) disk models end in gravitational runaway, while high-$`Q`$ ($`Q=1.5`$) disk models remain extremely stable. In both high and low $`Q`$ limits, the ultimate outcome is decided independent of the presence or absence of magnetic fields. When $`Q`$ is marginal ($`Q=1`$), on the other hand, we find that only magnetized models form bound condensations, suggesting that Parker instability may play a supplementary role in destabilization. The Parker instability raises the level of density fluctuation and may help channel power in the saturated state into small wavenumber modes. It is the “rejuvenated” swing amplification (see Paper I) that ultimately drives the system to form dense condensations, however. For cases in which bound clouds do form, the sizes and spacings are not characteristic of Parker modes, but comparable to the results from thin disk models (with typical masses $`M(23)`$ times the thin disk Jeans mass).

4. Weakly shearing thick disks, like thin disks, are susceptible to the magneto-Jeans instability (MJI), in which tension forces from azimuthal magnetic fields vitiate the stabilizing Coriolis force. The azimuthal wavelength criterion for MJI in finite-thickness disks is similar to that of the Parker instability ($`k_yH\stackrel{<}{}1`$), and modal growth rates are also similar. Our nonlinear simulations show, however, that MJI under weak-shear (spiral arm) conditions is ultimately much more violent than the Parker instability, and can form dense condensations within one orbital time. This is because the Parker instability is self-limiting, whereas the MJI is a runaway process. The mass of dense condensations amounts typically to the thick-disk Jeans mass ($``$ 4 times the thin-disk Jeans mass).

### 5.2 Discussion

The numerical simulations in this paper suggest that formation of bound clouds cannot occur as a direct consequence of Parker instability, although Parker instability may play a role in seeding swing amplification in high-shear environments or MJI in regions of weak shear. There are several reasons why the Parker instability is ineffective at producing large, dense clouds. First, although the growth time of the Parker instability from linear-theory analyses is a few times $`10^7`$ yrs (cf, Kim & Hong 1998), the process does not proceed in a runaway fashion, but is stabilized at a relatively modest density by magnetic tension forces (e.g., Mouschovias 1974). Second, while the undular Parker instability prefers wavelengths in the azimuthal direction $`\lambda _y2\pi H`$ (i.e., $`k_yH1`$ in rough agreement with observed spacing between molecular complexes), there is no similar preferred scale in the radial direction (e.g., Parker 1967; Giz & Shu 1993; J. Kim et al. 1998). Growth rates increase slowly with $`k_x`$, and are nearly uniform for $`k_xH\stackrel{>}{}5`$. Indeed, even though the initial perturbation amplitudes in our simulations increase toward large scale, we find that small-$`\lambda _x`$ modes dominate the evolution. The results is long filamentary structures with a radial spacing about 10 times smaller than the azimuthal spacing in surface density maps during the initial growing phase (see Fig. 6).121212The corresponding radial wavelengths are one or two times the disk scale height for symmetric and antisymmetric structures, respectively, possibly indicating nonlinear selection of these modes. Because the perturbed self-gravitational potential $`\delta \mathrm{\Phi }_s\delta \rho /k^2`$ is increasingly small for large $`k_x`$, effects of self-gravity on development of the Parker-Jeans instability are weak for the dominant nonlinear modes.

Third, while two-dimensional simulations for the Parker instability inarguably reach Mouschovias (1974)-type final static equilibrium configurations (cf, Santillán et al. 2000), three-dimensional simulations (J. Kim et al. 1998, 2000, and this work) indicate that such equilibrium configurations are transient at best, evolving rapidly into a state where density is relatively smooth in the horizontal direction. It is unclear whether this transition is initiated by an instability of the static undulating configurations, as claimed by Lachièze-Rey et al. (1980), by active magnetic reconnection near the midplane caused by small-scale chaotic motions, or simply by comparable nonlinear saturation levels of many high-$`k_x`$ modes. When shear is present, kinematically driven phase mixing contributes to washing out large-scale radial structures. In any case, the three-dimensional simulations show that it is almost impossible to find any signature of the undular Parker instability in the later stages of evolution.

The result that Parker instability cannot directly account for GMC formation has previously been emphasized by other authors. Kim, Ryu, & Jones (2001), for example, find only a maximum factor of 1.5 enhancement in surface density in their three-dimensional non-self-gravitating simulations with solid-body rotation, and conclude that Parker mode alone cannot be responsible for forming molecular clouds (see also e.g., Basu, Mouschovias, & Paleologou 1997; J. Kim et al. 1998, 2000; Santillán et al. 2000). These authors speculate on the possibility of the Parker instability coupled with other physical processes such as self-gravity, differential rotation, radiative cooling, and spiral arm potentials. We show in the present paper that self-gravity and galactic differential rotation are not much help. Therefore, we argue that the Parker instability, even combined with gaseous self-gravity and differential rotation, cannot be regarded as the primary formation mechanism for giant molecular clouds at least in galactic disks at large. While effects on the Parker instability of stronger self-gravity and rapidly varying shear inside spiral arms have yet to be studied in detail, we do not anticipate significant differences in the main conclusions. Although the width of a spiral arm might help in promoting Parker modes of similar $`\lambda _x`$ over shorter-wavelength competitors (cf., Franco et al. 2002), it is not yet clear how effective this can be if the arm transit time is comparable to the growth time.

While larger or smaller shear rates do not seem to alter the basic character of Parker instabilities appreciably, the shear rate is crucial for selecting which type of (in-plane) gravitating amplifier can operate. In a weakly shearing regions, swing is absent, but the MJI is very powerful at initiating self-gravitational condensations. In Paper I, we showed MJI is a potential means to prompt observed active star formation near the central parts of galaxies where rotation curves are nearly rigid-body (see also Elmegreen 1987). In Paper II, we argued that observed spur structures jutting radially outwards from spiral arms can be understood as a direct consequence of the MJI occurring within the arms. A spiral arm potential has three important effects on the gaseous medium. (1) By compressing gas, the Jeans length – which determines the preferred MJI scale along the mean magnetic field – is reduced; the growth time for self-gravitating modes also drops. (2) The spatially-varying density distribution causes a gradient in the shear profile perpendicular to the arm; since transiting gas first experiences reversed shear and then returns to normal shear, the spatially-varying profile tends to keep radial wavenumbers small within the arm, which enhances the efficiency of MJI. (3) It reduces the radial scale length for mass collection by MJI to at most the arm width, so that the masses of clumps produced are limited.

Compared to our previous results from razor-thin disks, the weaker self-gravity in three dimensions could imply somewhat more massive ($`10^7\mathrm{M}_{}`$) clouds would tend to form from MJI under realistic conditions. On the other hand, the weaker self-gravity may allow stronger spiral potential perturbations that give a higher density contrast between arm and interarm regions in stable quasi-axisymmetric spiral shock profiles (cf. Paper II). The Jeans mass at the spiral arm density peak could be correspondingly smaller, producing lower mass clouds in three-dimensional spiral arm regions. Since the marginal wavenumber of MJI parallel to $`𝐁`$ is about the same as that of the Parker instability in spiral arm conditions (see §4.3), one might expect some cooperation between MJI and the Parker instability inside spiral arms. The present models under uniform in-plane background conditions do not point to significant modifications of the nonlinear MJI by Parker modes. Effects from Parker-MJI interactions could be more prominent with a realistic spiral arm model (cf. Martos & Cox 1998), however, and it would be particularly interesting to learn whether the final mass scale changes markedly compared to predictions from uniform disk models.

By focusing on cases with relatively strong ($`\beta `$=1) equilibrium azimuthal magnetic fields in our MHD models, we have concentrated solely on effects of the Parker instability and MJI in the present work. An additional dynamical process prevalent in magnetized disk systems is the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998). Because MRIs are most vigorous in their small-scale growth when the thermal pressure is large compared to magnetic pressure ($`\beta 1`$), they have been studied primarily in the context of accretion disks. However, Sellwood and Balbus (1999) have argued that the MRI may also be important in producing turbulence in galactic disks where the effective value of $`\beta `$ is close to unity, citing as evidence the significant nonthermal velocity dispersions ($`6\mathrm{km}\mathrm{s}^1`$) present in outer H I disks with low local star formation rates and apparent gravitational stability. It is of great interest to investigate whether coupling between the MRI and self-gravitating modes under galactic conditions could represent an important new mechanism for forming large-scale ISM condensations.

We are grateful to an anonymous referee and N. Turner for thoughtful comments and to S. White, T. Sano, and J. Lee for help in using IDL to make three-dimensional visualization. This work was supported by NASA grant NAG 59167.