format: text

text:
# Formation and Fragmentation of Gaseous Spurs in Spiral Galaxies

## 1 Introduction

Even the grandest of grand-design spiral galaxies abound with substructure overlying and interlaced with the dominant twin-armed global pattern. This substructure may take the form of feathering between the main arms, or of spurs or larger branches jutting nearly perpendicularly from the primary arms, then sweeping back in the same sense at larger radii (e.g., Lynds 1970; Weaver 1970; Elmegreen 1980; Roberts, Lowe, & Adler 1990). Secondary structures often appear in combination (e.g., Piddington 1973), and are sometimes associated with chains of interarm H II regions (van der Kruit & de Bruyn, 1976; Elmegreen, 1979). These intermediate-scale features occur in multiple-armed spirals as well as grand design types, but blend more unobtrusively into the complex overall structure. Even at quite small (circumnuclear) scales, special image processing techniques can reveal regular substructure in spiral arms (Lou et al, 2001). From morphological evidence in multiple colors combined with kinematic arguments, Elmegreen (1980) concluded that many spurs represent long-lived wavelike phenomena.

Although spurs have long been recognized as characteristic features from ground-based studies, the higher resolution afforded by space-based platforms opens possibilities for observing these structures in exquisite and unprecedented detail. In particular, the recently-released Hubble Heritage image of the Whirlpool Galaxy M51 (Scoville & Rector, 2001; Scoville et al, 2001) stunningly reveals that spurs project, at quasi-regular intervals, from essentially the whole length of the two main arms. These spurs are defined by narrow dust lanes, and are dotted with H II regions both near the main arms and well into the interarm regions.

What is the origin of all these feathers, spurs, and branches? Do they represent a stochastic, localized phenomenon, or does the global spiral pattern play a role in their formation? Small-scale, irregularly-oriented and -separated feathering could in principle arise in many ways. Proposals include gravitational induction of quasi-steady trailing responses to orbiting local mass condensations (Julian & Toomre, 1966; Byrd, 1983; Byrd, Smith, & Smith, 1984), transient swing amplification of “shearing bits and pieces” from a near-uniform interarm background (Goldreich & Lynden-Bell, 1965; Julian & Toomre, 1966; Toomre, 1981), shocking of the interstellar medium (ISM) produced by ballistic clouds traversing spiral arms (Pikelner, 1970; Kaplan & Pikel’ner, 1974), and percolation in a shearing environment of self-propagating star formation (e.g., Gerola & Seiden 1978; Schulman & Seiden 1986). Large-scale interarm branches that lie in relative radial isolation could be associated with ultraharmonic resonances of the primary pattern’s potential (Shu, Milione, & Roberts, 1973; Elmegreen, Elmegreen, & Seiden, 1989; Elmegreen & Elmegreen, 1990).

In addition to structural features for which the large-scale spiral pattern is irrelevant, and those for which it is crucial (via a resonance), there may be features that are organized in a spatially-coherent plan by the spiral arms, yet which develop similarly to local, self-gravitating, shearing wavelets. We believe that the spurs recently observed in M51 (and likely to show themselves in other galaxies upon high-resolution examination) represent this sort of hybrid phenomenon.

The idea that spiral arms may trigger the growth of smaller-scale, self-gravitating structures – particularly in the gaseous medium – is long established, and supported by both observational and theoretical studies. Observational evidence includes the consistency with characteristic Jeans-unstable values of observed masses and separations of H I superclouds (Elmegreen & Elmegreen, 1983; Elmegreen, 1987b; Knapen et al, 1993) and giant molecular associations (Vogel et al, 1988; Rand, 1993; Sakamoto, 1996; Thornley & Mundy, 1997a, b; Sakamoto et al, 1999). In addition, H II regions and star complexes are often distributed along arms in a “beads on a string” pattern, with a separation of roughly 3 times ($`14`$ kpc) the full arm thickness, independent of pitch angle and galactocentric radius (Elmegreen & Elmegreen, 1983). Elmegreen (1994) argued that this spacing of giant H II regions is attributable to gravitational collapse of magnetized gas within spiral arms.

An origin of spur formation with gravitational instability in spiral arms is natural because it relies only upon the fundamental physical agents known to be present: galactic differential rotation, compression by the spiral pattern, and self-gravity; we shall show that well-developed spurs may require magnetic fields as well. Two barriers to gas condensations on a large scale are insufficient surface density and excessive shear; because both of these impediments are reduced simultaneously in spiral arms (e.g., Roberts 1969; Elmegreen 1987a), these regions are prone to self-gravitating instabilities that can produce nonaxisymmetric structure. Condensations may also grow inside spiral arms under the combined action of the Parker and gravitational instabilities, subsequently being carried downstream by galactic differential rotation to appear as spurs. In this paper, however, we isolate and concentrate only on local gravitational instability occurring in two-dimensional (2D), razor-thin disks; the effects of the Parker instability will be studied in future work.

Several previous studies have considered self-gravitating growth of perturbations within spiral arms in the linear regime. Balbus & Cowie (1985) analyzed local stability of quasi-axisymmetric disturbances (i.e. $`𝐤`$ perpendicular to the local spiral arm), and showed that in the linear regime growth is intrinsically transient, stabilized by the expanding, shearing flow as the medium leaves the arm region. Balbus (1988, hereafter B88) analyzed local amplification of shearing wavelets with arbitrary initial local wavenumber, allowing for a varying background shear profile self-consistent with the varying ISM compression in the arm/interarm regions. Elmegreen (1994) incorporated effects of magnetic fields in a WKB analysis of instantaneous growth in a spiral arm region, without explicitly taking into account the expansion and varying shear flow downstream (relevant to the swing mechanism; see B88).

Of previous work, the most detailed study relevant to the problem of spiral arm spur formation is that of B88. He showed that transient gravitational instabilities within compressed regions have two preferred directions for initial wavenumbers: nearly parallel and nearly perpendicular to the arm. He also suggested that closely-spaced periodic spurs may naturally develop from the former, while the latter could produce (non-periodic) spurs under certain conditions for the background velocity field.

In this paper, we investigate nonlinear evolution of local self-gravitating perturbations in the magnetized, gaseous ISM as they interact with the spiral potential in a differentially-rotating galaxy. Our study extends the work of B88 by including the effects of magnetic fields and nonlinearity. It extends our own previous nonlinear study of self-gravitating evolution in shearing, magnetized, featureless disks (Kim & Ostriker, 2001) to incorporate the strong density and velocity gradients in the equilibrium gas profile driven by a stellar spiral arm potential. For thin, featureless disks, there are two different mechanisms that grow density structure: the magneto-Jeans instability when shear is weak or magnetic fields are strong; and the magnetically-modified swing amplifier when shear is strong and the magnetic field is moderate or weak. Here, we study the modifications of these mechanisms by the presence of the spiral arm.

We proceed technically by considering a local patch of an infinitesimally thin gaseous disk subject to an external (stellar plus dark matter) potential. The model disk is assumed to be isothermal and initially magnetized parallel to the spiral arms. We adopt Roberts (1969)’s coordinate system having two orthogonal axes parallel and perpendicular to the gradient in the local external spiral potential. We first obtain one-dimensional (1D) density and velocity profiles with variation only in the direction perpendicular to the arm. We then apply small perturbations and follow their development with 2D direct numerical simulation. For comparison with these simulations, we also integrate the magnetized equivalent of the linear shearing-wavelet equations of B88.

The remainder of this paper is organized as follows: In §2, we introduce the Roberts coordinate system, present the vertically integrated, 2D magnetohydrodynamic (MHD) equations in these spiral-arm coordinates, and describe our simulation model parameters and numerical method. In §3, we pursue 1D spiral shock equilibria under the given background conditions, and study quasi-axisymmetric instabilities. We present in §4 a set of MHD simulations with zero, sub-equipartition, and equipartition magnetic fields, to explore 2D structure formation. We analyze density and velocity structures of the spurs that we find form, and compare these structures with analytic predictions. We present density snapshots from late stages of evolution for several models, showing that spurs frequently fragment to produce dense knots with masses comparable to giant clouds or massive young stellar associations. We summarize our results and discuss the conclusion of present work in §5. In Appendix A, we present the linear-theory perturbation equations for magnetized flows in spiral arms.

## 2 Basic Equations and Model Parameters

### 2.1 Coordinates and Basic Equations

We study the nonlinear response of a self-gravitating, magnetized, differentially rotating, razor-thin, gaseous disk to an externally-imposed spiral potential. The spiral potential, which is thought of as arising primarily from the stellar component, is assumed to be tightly wound with a pitch angle $`i=\mathrm{arctan}(\lambda _R/\lambda _\theta )1`$ and rigidly rotating relative to an inertial frame at a constant pattern speed $`\mathrm{\Omega }_p`$. The spiral forcing tends to induce noncircular motion in gaseous disks that would otherwise rotate with a angular velocity $`\mathrm{\Omega }(R)`$ in the azimuthal direction $`\widehat{\theta }`$, where $`R`$ is the galactocentric radius. The induced radial velocities are generally only a few percent of $`R\mathrm{\Omega }(R)`$, but enough to compress gas into a highly nonlinear state. We investigate the gravitational instabilities that may develop in these highly-compressed spiral-arm regions.

Following Roberts (1969) and B88, we construct a local Cartesian frame whose center lies at a galactocentric distance $`R_0`$ and rotates at $`\mathrm{\Omega }_p`$ about the galactic center. The two coordinates correspond to radial $`(RR_0)`$ and angular ($`R_0(\theta \mathrm{\Omega }_pt)`$) displacements with respect to the center of our moving frame. We then tilt our frame by an angle $`i`$ such that $`\widehat{x}`$ and $`\widehat{y}`$ correspond to the directions perpendicular and parallel to the local spiral arm, respectively, and consider a local rectangular box with size $`L_x\times L_y`$, as depicted in Figure 1. We expand the compressible, ideal MHD equations in the rotating frame, and take the approximations that $`\mathrm{sin}i1`$ (tightly wound arms), $`x,yR_0`$ (local model), and the velocities induced by the spiral arm are much smaller than $`R_0\mathrm{\Omega }(R_0)`$. Neglecting the terms arising from the curvilinear geometry and integrating the resulting equations in the vertical direction, we obtain the following set of 2D equations:

$$\frac{\mathrm{\Sigma }}{t}+(\mathrm{\Sigma }𝐯_T)=0,$$
(1)
$`{\displaystyle \frac{𝐯}{t}}+𝐯_T𝐯`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Sigma }}}\mathrm{\Pi }+{\displaystyle \frac{1}{4\pi \mathrm{\Sigma }}}(\times 𝐁)\times 𝐁`$ ()
$`+q_0\mathrm{\Omega }_0v_x\widehat{𝐲}2𝛀_0\times 𝐯(\mathrm{\Phi }_g+\mathrm{\Phi }_{\mathrm{ext}}),`$
$$\frac{𝐁}{t}=\times (𝐯_T\times 𝐁),$$
(3)
$$^2\mathrm{\Phi }_g=4\pi G\delta (z)\mathrm{\Sigma },$$
(4)
and

$$\mathrm{\Pi }=c_s^2\mathrm{\Sigma },$$
(5)
(cf, Roberts (1969); Roberts & Yuan (1970); Shu, Milione, & Roberts (1973); Balbus & Cowie (1985); B88). Here, $`\mathrm{\Sigma }`$ is the surface density, $`\mathrm{\Pi }`$ is the 2D vertically-integrated pressure, $`\mathrm{\Phi }_g`$ is the self-gravitational potential, and $`𝐁`$ is the midplane value of the 3D magnetic field times the square root of the unperturbed ratio of surface density to midplane volume density. The total velocity as viewed in the rotating frame is denoted by $`𝐯_T`$, while $`𝐯`$ represents the part induced by the external spiral potential $`\mathrm{\Phi }_{\mathrm{ext}}`$; that is, $`𝐯𝐯_T𝐯_0`$, where $`𝐯_0`$ is the velocity arising from galactic rotation111By direct coordinate transformation, $`𝐯_0=[R_0(\mathrm{\Omega }_0\mathrm{\Omega }_p)q_0\mathrm{\Omega }_0(x\mathrm{cos}iy\mathrm{sin}i)](\mathrm{sin}i\widehat{𝐱}+\mathrm{cos}i\widehat{𝐲})`$. Within the local, tightly-winding limit ($`x,yR_0`$ and $`\mathrm{sin}i1`$) that we adopt, the background circular velocity reduces to $`𝐯_0=R_0(\mathrm{\Omega }_0\mathrm{\Omega }_p)\mathrm{sin}i\widehat{𝐱}+[R_0(\mathrm{\Omega }_0\mathrm{\Omega }_p)q_0\mathrm{\Omega }_0x]\widehat{𝐲}`$. This form of the background velocity permits local shearing periodic boundary conditions.. In equation (2), $`\mathrm{\Omega }_0\mathrm{\Omega }(R_0)`$ and $`q_0(d\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}R)|_{R_0}`$ measures shear rate in the background flow in the absence of a spiral pattern, while $`c_s`$ in equation (5) is the isothermal sound speed. In terms of $`q_0`$ and $`\mathrm{\Omega }_0`$, the local epicyclic frequency $`\kappa _0`$ is given by $`\kappa _0^2R^3d(R^4\mathrm{\Omega }^2)/dR|_{R_0}=(42q_0)\mathrm{\Omega }_0^2`$. Finally, $`G`$ and $`\delta `$ in equation (4) are the gravitational constant and the Kronecker delta, respectively.

As expressed by equation (5), in this paper we adopt an isothermal equation of state. In equations (2) and (3), we treat the effective scale height of the magnetic field distribution as a constant in both space and time. In the present height-integrated formulation, our dynamical model does not capture the potential consequences of Parker instability and magnetorotational instabilities for disk evolution (cf. Kim & Ostriker 2001). Nevertheless, because our equations contain the essential physical ingredients involved in self-gravitational instabilities inside spiral arms, the present work represents the first step in a more comprehensive three-dimensional study.

To complete the governing equations, we must specify the external spiral potential $`\mathrm{\Phi }_{\mathrm{ext}}`$. For compatibility with our local model, $`\mathrm{\Phi }_{\mathrm{ext}}`$ must be periodic in $`x`$ and independent of $`y`$; we adopt the following simple form:

$$\mathrm{\Phi }_{\mathrm{ext}}=\mathrm{\Phi }_0\mathrm{cos}\left(\frac{2\pi x}{L_x}\right),$$
(6)
which is a local analog of a logarithmic potential of Roberts (1969) and Shu, Milione, & Roberts (1973). Since $`x`$ varies from $`L_x/2`$ to $`L_x/2`$, $`\mathrm{\Phi }_{\mathrm{ext}}`$ attains its minimum at the center ($`x=0`$) for $`\mathrm{\Phi }_0<0`$. The distance $`L_x`$ is the arm-to-arm separation, equal to $`2\pi R_0\mathrm{sin}i/m`$ for an $`m`$-armed spiral.

### 2.2 Model Parameters

In this paper, we adopt an isothermal equation of state; for scaling our solutions, we shall use an effective isothermal speed of sound $`c_s=7\mathrm{km}\mathrm{s}^1`$, corresponding to a mean 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). For the spiral potential parameterization, we take pattern speed $`\mathrm{\Omega }_p=\mathrm{\Omega }_0/2`$ and pitch angle $`\mathrm{sin}i=0.1`$, respectively. For a two-armed spiral pattern located at distance $`R_0=10`$ kpc from the galactic center, we thus have $`L_x=\pi R_0\mathrm{sin}i=3.1`$ kpc. In the solar neighborhood, the angular velocity of the Galactic rotation is $`\mathrm{\Omega }_0=26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1`$ (Binney & Tremaine, 1987), and $`q_01`$ for a near-flat rotation, so that $`\kappa _02^{1/2}\mathrm{\Omega }_0`$. The corresponding orbital period is $`t_{\mathrm{orb}}2\pi /\mathrm{\Omega }_0=2.4\times 10^8\mathrm{yr}(\mathrm{\Omega }_0/26\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1)^1`$. The size of the simulation domain in the $`y`$-direction is chosen to be $`L_y=2L_x`$ (we find that our results are generally independent of this choice).

In the absence of external potential perturbation, our model disks have uniform surface density $`\mathrm{\Sigma }_0`$, and uniform magnetic field $`𝐁_0`$ that points in the $`\widehat{y}`$-direction. To characterize $`\mathrm{\Sigma }_0`$ and $`B_0`$, we introduce two dimensionless parameters,

$$Q_0\frac{\kappa _0c_s}{\pi G\mathrm{\Sigma }_0}$$
(7)
and

$$\beta _0\frac{c_s^2}{v_\mathrm{A}^2}=\frac{4\pi \rho _0c_s^2}{B_{0,3\mathrm{D}}^2},$$
(8)
where the Alfvén speed is defined by $`v_\mathrm{A}B_0^2/4\pi \mathrm{\Sigma }_0B_{0,3\mathrm{D}}^2/4\pi \rho _0`$ with 3D field strength $`B_{0,3\mathrm{D}}`$ at the disk midplane. The Toomre $`Q_0`$ parameter is a measure of the axisymmetric stability of a rotating, unmagnetized disk of uniform density. For the background conditions adopted in this paper, the corresponding mean gas surface density and mean magnetic field strength are

$$\mathrm{\Sigma }_0=\frac{19}{Q_0}\mathrm{M}_{}\mathrm{pc}^2\left(\frac{c_s}{7.0\mathrm{km}\mathrm{s}^1}\right)\left(\frac{\kappa _0}{36\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1}\right),$$

and

$$B_{0,3\mathrm{D}}=\frac{3.2}{\sqrt{\beta _0}}\mu \mathrm{G}\left(\frac{n_\mathrm{H}}{1\mathrm{cm}^3}\right)^{1/2}\left(\frac{c_s}{7.0\mathrm{km}\mathrm{s}^1}\right).$$

Note that the total mass contained in the box is then roughly $`M_{\mathrm{tot}}=L_xL_y\mathrm{\Sigma }_0(4/Q_0)\times 10^8\mathrm{M}_{}.`$

Finally, the strength of the external spiral potential is characterized by

$$F\frac{2}{\mathrm{sin}i}\left(\frac{|\mathrm{\Phi }_0|}{\mathrm{\Omega }_0^2R_0^2}\right),$$
(9)
which measures the amplitude of the perturbed radial force $`2\pi |\mathrm{\Phi }_0|/L_x`$ as a fraction of the mean axisymmetric gravitational force field (Roberts, 1969). We vary $`Q_0`$, $`\beta _0`$, and $`F`$ to represent disks in various physical conditions, while fixing the other input parameters as described above. Note that other choices of $`\mathrm{\Omega }_0`$, $`c_s`$, $`L_x`$, and $`\mathrm{sin}i`$ may be substituted for our fiducial values, provided that the ratio $`c_s\mathrm{\Omega }_0/L_x`$ remains unchanged, and $`F`$ is interpreted as varying proportional to $`\mathrm{sin}i`$.

### 2.3 Numerical Method

The nonlinear evolution of gaseous model disks is followed by integrating the governing set of equations (1)$``$(5) 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. 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 this work we implement shearing box boundary conditions, in which the $`y`$-boundaries are perfectly periodic and the $`x`$-boundaries are shearing-periodic (Hawley, Gammie, & Balbus, 1995). For gaseous self-gravity, we implement a shearing-sheet Poisson solver (Gammie, 2001) via FFTs such that $`\mathrm{\Phi }_k=2\pi G\mathrm{\Sigma }_k/|k|.`$ For less diffusive advection of hydrodynamic variables, we apply the velocity decomposition method which treats the contribution from the background shearing parts as source terms. 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. Our implementation of the ZEUS code has been fully tested on a variety of test problems and used for the study of gravitational instabilities in galactic disks without spiral-arm features (Kim & Ostriker, 2001). For the interested reader, we refer to Kim & Ostriker (2001) for a detailed description of the code and its test results. Our resolution for the 2D simulations presented in this paper is 256$`\times `$512 zones.

## 3 One-Dimensional Spiral Shocks and Their Stability

It is straightforward to use our numerical code to obtain 1D solutions for the gaseous response to an applied spiral potential. In this section, we present radially-localized solutions for steady-state 1D spiral shock structures and investigate their stability by performing fully nonlinear, time-dependent calculations. Steady-state 1D solutions including shocks were originally obtained by Roberts (1969), Roberts & Yuan (1970), Shu et al (1972), Shu, Milione, & Roberts (1973), and Tomisaka (1987) for non-self-gravitating flows, and by Lubow, Balbus, & Cowie (1986) for viscous, self-gravitating, unmagnetized flows. Woodward (1975) followed a time-dependent approach and showed that spiral shocks indeed form within one or two crossing times of gas through the spiral density waves. While his treatment allowed for the effects of ultraharmonic resonances (cf. Shu, Milione, & Roberts 1973) in driving the nonlinear gaseous response, Woodward (1975) did not address the stability of spiral shocks when self-gravity and magnetic fields are incorporated. Similar works for 2D galactic shocks including both spiral and bar potentials were reported by Roberts, Huntley, & van Albada (1979) and van Albada & Roberts (1981).

Linear analyses of 1D gravitational instabilities of the ISM inside spiral arms were performed by Balbus & Cowie (1985); here we shall present the results of nonlinear simulations. Balbus & Cowie (1985) focused on the development histories of gravitationally-amplified perturbations, showing that in linear theory the background expansion of the flow as it exits the arm limits growth. Thus, in analogy to the situation for nonaxisymmetric perturbations in a shearing background (Goldreich & Lynden-Bell, 1965), there are no true linear instabilities but only a transient growing phase, analogous to the swing amplifier. Just as nonlinear effects nevertheless lead to an instability threshold for the swing amplifier (Kim & Ostriker, 2001), one may expect a nonlinear instability threshold in the 1D spiral shock problem. In this section, we concentrate on finding critical $`Q_0`$ values for the existence of equilibrium spiral shocks for given galactic conditions. As we vary only the magnetic field strength and the spiral perturbation amplitudes, our coverage of the parameter space is of course incomplete. The trends shown in the critical $`Q_0`$ curves are interesting in themselves, but are also important for delimiting the portion of parameter space for which nonaxisymmetric studies are warranted. One-dimensional equilibrium shock profiles constructed in this section will serve as initial conditions for our nonaxisymmetric simulations in §4.

We begin by considering a non-self-gravitating medium with uniform density $`\mathrm{\Sigma }_0`$, background shear profile $`𝐯_0`$, and uniform magnetic field characterized by $`\beta _0`$. We then impose an external spiral perturbation and slowly increase its amplitude up to a desired level, $`F`$ (eq. ). The spiral forcing drives magnetosonic waves, but these waves with $`c_s/v_{0,x}<1`$ soon steepen, eventually developing a shock front. The surface density first overshoots and oscillates about the corresponding steady-state solution for given $`\beta _0`$ and $`F`$, and then gradually converges to it. The typical time to reach a stationary state is about 5 orbital times. We check the code accuracy by comparing the numerical solutions with published solutions obtained from solving, with a shock-fitting procedure, the time-independent ordinary differential equations (Roberts, 1969; Shu et al, 1972; Shu, Milione, & Roberts, 1973).

We then slowly turn on the self-gravitational force, as we did for $`F`$, to obtain a self-gravitating equilibrium profile corresponding to a given $`Q_0`$. Two examples of stable equilibria for parameter values $`Q_0=2.0`$, $`\beta _0=10`$, $`F=3\%`$ and $`Q_0=1.5`$, $`\beta _0=1,`$ $`F=3\%`$ are respectively illustrated in Figures 2 and 3 (solid curves). Note that gas is flowing from left to right. Dotted curves in Figures 2$`a`$-2$`e`$ and 3$`a`$-3$`e`$ indicate the unperturbed state with $`F=0`$, while the dashed lines in Figures 2$`b`$ and 3$`b`$ mark the sound speed. For comparison, we also plot the surface density profiles for the non-self-gravitating counterparts as dashed lines in Figures 2$`a`$ and 3$`a`$ . A well-defined shock front occurs at $`x/L_x0.02(0.05)`$ with a thickness $`\mathrm{\Delta }x/L_x0.03(0.02)`$ for the $`\beta _0=10(1)`$ case. Across the shock front, $`v_x`$ experiences a sharp deceleration from supersonic to subsonic magnitude, while the gas and magnetic field are both compressed. After the density peak, gas begins to expand and $`v_x`$ becomes again supersonic after crossing another sonic point at $`x/L_x0.079(0.086)`$ for $`\beta _0=10(1)`$. For flows with one independent spatial variable, mass and magnetic flux conservation guarantees $`\mathrm{\Sigma }Bv_x^1`$, so that $`\beta v_x`$, as confirmed by Figures 2$`b`$,2$`d`$ and 3$`b`$,3$`d`$. The conservation of the potential vorticity222 In these expressions for $`\xi `$, the curl operator and velocities are defined in the local Cartesian frame built in §2.1. When expressed in standard cylindrical coordinates retaining curvilinear terms, $`\xi =|\times 𝐕|/\mathrm{\Sigma }=|\times 𝐯_T+2𝛀_p|/\mathrm{\Sigma }`$, where $`𝐕`$ is the fluid velocity seen in the inertial frame and $`𝐯_T𝐕R\mathrm{\Omega }_p\widehat{\varphi }`$. $`\xi |\times 𝐯_T+2𝛀_0|/\mathrm{\Sigma }=|\times 𝐯+(2q_0)𝛀_0|/\mathrm{\Sigma }`$ (Hunter, 1964; Gammie, 1996, 2001) leads to an equilibrium profile $`Q=Q_0(\mathrm{\Sigma }/\mathrm{\Sigma }_0)^{0.5}`$ for isothermal gas or $`Q=Q_0\gamma ^{1/2}(\mathrm{\Sigma }/\mathrm{\Sigma }_0)^{\gamma /21}`$ for polytropic law $`\mathrm{\Pi }\mathrm{\Sigma }^\gamma `$ (Balbus & Cowie, 1985), and implies that the local shear rate varies as

$$q\frac{\mathrm{ln}\mathrm{\Omega }}{\mathrm{ln}R}|_{R_0}=\frac{1}{\mathrm{\Omega }_0}\frac{dv_{y,T}}{dx}=2(2q_0)\frac{\mathrm{\Sigma }}{\mathrm{\Sigma }_0},$$
(10)
(Kim & Ostriker, 2001). As a consequence, shear is reversed where $`\mathrm{\Sigma }/\mathrm{\Sigma }_0>2`$ for $`q_0=1`$, as shown in Figures 2$`c`$ and 3$`c`$.

Our simulation domain contains an ultraharmonic resonance where $`v_{x,0}^2c_s^2=(L_x/2\pi )^2\kappa ^2/n^2`$ for a non-self-gravitating, unmagnetized medium, with integer $`n`$ denoting the order of the resonance (Shu, Milione, & Roberts, 1973). Using the parameters adopted in §2.2, one can show that $`x/L_x=0.245`$ corresponds to a second harmonic resonance ($`n=2`$). Although this resonance is easily detectable when spiral arm forcing is relatively weak (Shu, Milione, & Roberts, 1973), the local maximum around $`x/L_x0.25`$ in the dashed curve for the non-self-gravitating model in Figure 2$`a`$ still traces the (weak) second harmonic resonance. Embedded (strong) magnetic fields and self-gravity significantly change the resonant conditions (cf. Lubow, Balbus, & Cowie 1986), so that $`\beta _0=1`$ models (e.g., Fig. 3$`a`$) and strongly self-gravitating $`\beta _0=10`$ models (e.g., solid line in Fig. 2$`a`$) do not show any indication of the resonance.

In general, self-gravity tends to move the shock front downstream, enhance the maximum surface density, and widen the shock front333Since the gravitational potential is continuous across the shock front, self-gravity does not affect the shock jump conditions (cf. Shu 1992). The peak density is not always reached immediately behind the shock (Shu, Milione, & Roberts, 1973). However, the equilibrium profiles presented in Figures 2 and 3 could not resolve the separation between the shock front and the locus of the maximum density.. Lubow, Balbus, & Cowie (1986) found that with viscosity included (so as to represent a nonzero mean free path in a cloud-fluid), the gas self-gravity suppresses the tendency of the gas to form a shock when the gas content is large. We find the similar result that when the spiral arm forcing is relatively weak (i.e., the gaseous self-gravity exceeds the force due to stars), the shock disappears, leaving quite symmetric density configurations (see, for example, model H1 in Fig. 5$`a`$). Notice the magnitudes of $`\mathrm{\Phi }_{\mathrm{ext}}`$ and $`\mathrm{\Phi }_g`$ are comparable in Figures 2$`f`$ and 3$`f`$, although $`\mathrm{\Phi }_g`$ has larger gradients.

With smaller $`Q_0`$ and/or larger $`F`$ (i.e., larger self and/or external gravity), spiral shock fronts move further downstream, and the arm-to-interarm contrast in the stationary density profile grows, increasing the susceptibility to gravitational instability. For models with $`Q_0`$ less than some critical value (for other parameters fixed), no stationary configuration can be found. This is either because the equilibrium becomes nonlinearly unstable for sufficiently strong self-gravity, or because no equilibrium for a given background condition is compatible with the chosen pattern speed and strength of the spiral potential. These critical $`Q_0`$ values together with corresponding local values of $`Q`$ at the density peak, $`Q_{\mathrm{sp}}Q_0(\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_0)^{1/2}`$, are indicated in Figure 4.

Gravitationally unstable cases, as represented by filled symbols in Figure 4, occur when the magnetic field and external forcing are weak. Stronger magnetic fields provide pressure support against gravitational collapse if $`Q_0`$ is not as small. As $`Q_0`$ is lowered, the equilibrium shock front moves gradually downstream, where the gradient of spiral potential is smaller. Beyond a certain point (when the shock approaches the maximum of $`\mathrm{\Phi }_{\mathrm{ext}}`$), the shock front is no longer stationary but moves back and forth, with rapidly changing density. This implies that steady-state solutions of equations (1)$``$(5) do not exist for a given background state. The boundaries of cases with this sort of behavior are marked by open symbols in Figure 4. Note that the critical $`Q_0`$ values, $`Q_{0,c}`$, for purely hydrodynamic models all represent instances of gravitational instabilities, ranging over $`Q_{0,c}1.02.8`$ for $`F6\%`$. As $`\beta _0`$ decreases (stronger mean magnetic fields), the range of $`Q_{0,c}`$ also decreases, giving $`Q_{0,c}0.92.0`$ for $`\beta _0=10`$ and $`Q_{0,c}0.70.9`$ for $`\beta _0=1`$.

In cases that are gravitationally unstable, perturbations grow as they move downstream off the arm. For this gravitational runaway to occur, amplification of any perturbations must be high enough to drive the system into the nonlinear regime before the onset of stabilization associated with physical expansion and high interarm $`Q`$. Balbus & Cowie (1985) suggested that $`Q_{\mathrm{sp}}`$ is the parameter that characterizes the gravitational response of gas to applied perturbations. For marginally stable models with $`1\%F6\%`$, Figure 4 shows that $`Q_{\mathrm{sp}}0.730.98`$ for $`\beta _0=\mathrm{}`$, which is in good agreement with Balbus & Cowie (1985)’s finding that $`Q_{\mathrm{sp}}<1`$ is required (but not sufficient) for significant growth of perturbation. The main reason why $`Q_{\mathrm{sp}}<1`$ does not necessarily produce axisymmetric instability is of course that the background surface density is not uniform; at least a region of the shortest unstable wavelength must have $`Q<1`$. For cases where marginality represents incompatibility of equilibria with prescribed background conditions, $`Q_{\mathrm{sp}}`$ converges to 0.5 for $`\beta _0=10`$ and to 0.4 for $`\beta _0=1`$ as $`F`$ increases. Note that $`Q_{\mathrm{sp}}`$ does not include the stabilizing effect of magnetic fields; including magnetic fields increases the effective $`Q`$ by a factor $`(1+\mathrm{\Sigma }/(\beta _0\mathrm{\Sigma }_0))^{1/2}`$, which is generally larger than unity at the limiting value of $`Q_{\mathrm{sp}}`$.

## 4 Two-Dimensional Simulations

In the previous section, we discussed the nonlinear stability of the gas flow through spiral arms to perturbations in which all spatial gradients point in the direction transverse to the spiral arm ($`\widehat{x}`$-direction). These may be thought of as quasi-axisymmetric instabilities. We now generalize to study the development of perturbations in which gradients may point in any direction with respect to the spiral arm. Disturbances where the gradients lie primarily along the arm ($`\widehat{y}`$-direction) may be thought of as quasi-azimuthal (nonaxisymmetric) perturbations.

For the initial conditions of our 2D simulations, we use the equilibrium shock profiles obtained from 1D calculations in §3. We consider cases with arm-interarm contrast in surface density less than a factor 12. Initial perturbations are realized by a Gaussian random density field with flat power for $`1kL_x/2\pi 64`$ and zero power for $`64<kL_x/2\pi `$. For the amplitude of perturbations, we measure the standard deviation $`ϵ_0`$ of the initial density fluctuation in real space, and fix $`ϵ_0=1\%`$. These adopted density perturbations do not attempt to represent realistic interstellar perturbations: in real galaxies, perturbations may have power-law spectra in both velocity and density fields with significantly larger amplitudes. The chosen amplitude and spectrum, nevertheless, allow us to monitor evolutionary behavior in both linear and nonlinear regimes, and seed the most dominant nonaxisymmetric mode of the instability in the simplest possible way.

The model parameters and nonaxisymmetric simulation results we present in this section are summarized in Table 1. Column (1) labels each run. Columns (2) and (3) list the two basic parameters of our model disks: the Toomre stability parameter $`Q_0`$ and the plasma parameter $`\beta _0`$ characterizing the magnetic field strength (see eqs.  and ). The amplitude of spiral arm forcing is given in terms of $`F`$ (see eq. ) in column (4). The peak surface density $`\mathrm{\Sigma }_{\mathrm{max}}`$ and local $`Q_{\mathrm{sp}}`$ of the 1D, stationary solution are listed, respectively, in columns (5) and (6), while column (7) gives the arm width determined at $`\mathrm{\Sigma }=(\mathrm{\Sigma }_{\mathrm{max}}+\mathrm{\Sigma }_{\mathrm{min}})/2`$. In the cases where spurs form, their mean spacing ($`\lambda _y`$) measured along the arm is given in columns (8) and (9) in the units of the arm-to-arm distance $`L_x`$ and the arm width $`W`$, respectively. Finally, column (10) gives the corresponding normalized wavenumber $`K_{y,\mathrm{max}}\lambda _{\mathrm{J},\mathrm{sp}}/\lambda _y=2\pi (\lambda _yk_{\mathrm{J},\mathrm{sp}})^1`$ of spurs, where the local Jeans wavenumber $`k_{\mathrm{J},\mathrm{sp}}`$ at the density peak is defined by

$$k_{\mathrm{J},\mathrm{sp}}\frac{2\pi }{\lambda _{\mathrm{J},\mathrm{sp}}}=\frac{2\pi G\mathrm{\Sigma }_{\mathrm{max}}}{c_s^2}.$$
(11)
Note, for comparison, that magneto-Jeans instability in a uniform disk with embedded magnetic fields would give $`K_{y,\mathrm{max}}=\frac{1}{2}`$ for $`\beta _01`$ and $`K_{y,\mathrm{max}}\frac{1}{2}\frac{3}{4}`$ for $`\beta _0\stackrel{>}{}0.1`$, while swing amplification generally favors $`K_{y,\mathrm{max}}0.150.4`$ for $`\beta _0\stackrel{>}{}1`$ (Kim & Ostriker, 2001).

### 4.1 Purely Hydrodynamic Models

For purely hydrodynamical simulations, we select models that are very close to marginal stability for a given arm strength $`F`$. Note that the critical $`Q_0`$ values for $`\beta =\mathrm{}`$ are 1.67, 2.01, and 2.26, for $`F=1`$%, 2%, and 3%, respectively. In Figure 5, we display the initial surface density profiles and the evolutionary histories of maximum surface densities for model H1$``$H3.

Generally speaking, unmagnetized configurations that are quasi-axisymmetrically stable are found to be stable also to nonaxisymmetric perturbations, unless $`Q_0`$ and $`F`$ are quite small. In models H2 and H3, with a realistic arm-interarm contrast in surface density of $`\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_{\mathrm{min}}=6.2`$ and 11.5, respectively, the nonaxisymmetric growth of density is so mild that no gravitationally bound structure forms within the simulation interval. The two main mechanisms that can produce local growth of perturbations in unmagnetized, razor-thin disks are quasi-axisymmetric instability and swing amplification (see B88, who refers to the quasi-axisymmetric modes as Jeans instability). Since our chosen initial conditions for nonaxisymmetric simulations are already stable to 1D quasi-axisymmetric instability, one can only expect perturbations to grow as wave crests swing from leading to trailing configurations. Figures 5$`b`$ shows that swing amplification in models H2 and H3 is in fact very weak, producing only small-amplitude fluctuations in the density. This is because the high density compression in spiral arm regions produces reversed shear (see eq.  with $`q_0=1`$ and Figs. 2$`c`$, 3$`c`$) under which condition classical swing amplification is essentially shut off. Another viable mechanism within spiral arms is divergence swing amplification that allows growth of leading wavelets in expanding background flows (B88). It is, however, very difficult to supply leading perturbations to spiral arms in a natural way, because the perturbations carried in from interarm regions are preferentially trailing. In interarm regions where $`\mathrm{\Sigma }/\mathrm{\Sigma }_0<2`$ so background shear is in the same sense as epicyclic motion, on the other hand, $`Q`$ is high enough to suppress the classical swing mechanism. For example, model H3 has $`Q2.8`$ at $`x/L_x=0.3`$.

Because of its smaller $`Q_0`$ and flatter density distribution, model H1 experiences stronger swing amplification than models H2 and H3. With smaller $`F`$, self-gravity smears out the spiral shock and makes the equilibrium profile more symmetric. Lower $`\mathrm{\Sigma }_{\mathrm{max}}`$ implies the sense of shear is “normal” in the major portion of a disk. Still, as can be seen in Figure 5$`b`$, growth of perturbations occurs so slowly that it requires four successive passages through spiral arm regions to amplify them significantly. Had we begun the evolution with higher perturbation amplitudes, less time would be required to attain a fully nonlinear state. As a result of swing amplification, model H1 forms three (weak) spurs in the direction perpendicular to the spiral arm, with an average separation of $`2.1`$ kpc. When allowed to evolve further, they fragment and collapse.

The spurs that form in model H1 are shown in Figure 6, where we plot the density at $`t/t_{\mathrm{orb}}=4.7`$, and for comparison the density in the stable model H3 at $`t/t_{\mathrm{orb}}=4.0`$. Warped by the nonlinear background shear and expansion flow off the arm, the local wavefronts defining the spurs in model H1 are well described by the kinematic formula provided by B88. Namely, the local tangent to the wavefront is given by $`dy/dx=k_x/k_y=𝒯`$ with

$$𝒯\frac{1}{}\left[\frac{\kappa _0^2\mathrm{\Sigma }_{\mathrm{max}}}{2\mathrm{\Omega }_0^2\mathrm{\Sigma }_0}\tau 2_0^\tau 𝑑\tau \frac{k_x(0)}{k_y}\right],$$
(12)
where $`\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }=v_{x,T}/v_{x,T,\mathrm{min}}`$ is the local surface density expansion factor, $`\tau \mathrm{\Omega }_0_{x_{\mathrm{sp}}}^xv_{x,T}^1𝑑x`$ is a dimensionless elapsed time that is measured from the shock location (or density peak), $`x_{\mathrm{sp}}`$, $`k_y=2\pi /\lambda _y`$, where $`\lambda _y`$ is the spur spacing, and $`k_x(0)`$ is the local $`x`$-wavenumber at $`\tau =0`$. The function $`𝒯`$ is defined in terms of the unperturbed state. From mass conservation, we can write the time element in terms of the expansion factor $``$ as

$$d\tau =\frac{\mathrm{\Omega }_0\mathrm{\Sigma }_{\mathrm{max}}dx}{R_0(\mathrm{\Omega }_0\mathrm{\Omega }_p)\mathrm{sin}i\mathrm{\Sigma }_0},$$
(13)
so that $`𝒯`$ only depends on the variation of surface density with $`x`$ via $`(x)`$. Using equation (13), one can show that the second term in the square brackets in equation (12) is directly proportional to $`xx_{\mathrm{sp}}`$; the first term, from the definition of $`\kappa _0^2`$, equals $`(2q_0)\tau \mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_0`$. The curves whose local gradient is defined by $`𝒯`$ with $`k_x(0)/k_y=2.2`$ are overlaid on the spurs in Figure 6, showing excellent agreement with the results of numerical simulations.

Our numerical evidence of the relative stability of purely hydrodynamic systems to nonaxisymmetric perturbations might seem in conflict with B88’s conclusion from linear theory that nonaxisymmetric wavelets can grow rapidly within $`0.3`$ orbital times. In fact, our results are entirely consistent with the predictions of B88’s general framework, but our numerical models explore a different part of the parameter space of background conditions – especially $`Q_{\mathrm{sp}}`$ – from that explored by B88. In particular, B88 considered a very soft polytropic equation of state ($`\gamma =0.5`$), such that he could easily achieve $`Q_{\mathrm{sp}}=Q_0\gamma (\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_0)^{\gamma /21}0.6`$ from solar neighborhood conditions, while our isothermal ($`\gamma =1`$) models have $`Q_{\mathrm{sp}}>0.9`$ as Table 1 indicates. As Figure 5 in B88 shows, amplification factors are very sensitive to $`Q_{\mathrm{sp}}`$, dropping by about one order of magnitude as $`Q_{\mathrm{sp}}`$ increases from 0.52 to 0.67. Naive extrapolation of this result suggests that models with $`Q_{\mathrm{sp}}>0.9`$ may remain stable; we have indeed verified that the solutions of the corresponding linear equations (see Appendix A for magnetized versions) show very low amplifications for our background models. Because our 2D simulations require an initial equilibrium profile (rather than an arbitrary shape for the surface density, as is permissible for linear analyses), and because such 1D equilibria do not exist for very small $`Q_{\mathrm{sp}}`$, our 2D simulations do not cover the parts of parameter space in which $`Q_{\mathrm{sp}}`$ is small studied by B88.

Our hydrodynamic models are by no means exhaustive, and the uncertainties in defining a realistic background model leave open many possibilities. Nevertheless, both the models shown here and others we have studied support the conclusion that without magnetic effects, it is difficult to form large-scale spur-like structures in a gaseous arm that is not violently unstable to quasi-axisymmetric modes, and hence unlikely to represent a realistic equilibrium in the first place.

### 4.2 MHD Simulations with Sub-equipartition Magnetic Fields

For simulations with mean magnetic field strength below thermal equipartition, we select three $`\beta _0=10`$ (i.e., $`v_\mathrm{A}=c_s/\sqrt{10}`$) models that are stable to quasi-axisymmetric perturbations. The critical $`Q_0`$ values for $`\beta _0=10`$ are 1.58, 1.78, and 1.89, for $`F=1`$%, 2%, and 3%, respectively. Figure 7 shows the initial density distributions and the evolution of maximum surface densities for model MS1$``$MS3.

Compared to Figure 5, Figure 7 demonstrates magnetic destabilization of modes with $`𝐤`$ along the arm: perturbations in all three models grow very rapidly, leading to the formation of spurs that will eventually fragment and collapse. This destabilization in spiral arms may be understood in relationship to the nonaxisymmetric instabilities present in shearing, featureless disks. In addition to swing amplification, for which magnetic fields play a stabilizing role, magnetized thin disks are also subject to magneto-Jeans instability (MJI), in which tension forces from embedded magnetic fields resist the stabilizing Coriolis force (Elmegreen, 1987a, 1994; Kim & Ostriker, 2001). The efficiency of the MJI depends on the local field strength and the local shear rate. Even with weak magnetic fields, significant growth of perturbations can occur, provided that their azimuthal wavelength is less than the local Jeans wavelength and the local shear rate is not too large (see, e.g., Fig. 1 of Kim & Ostriker 2001). Although swing amplification may also contribute at a lower level, background spiral-arm conditions in models MS2 and MS3 provide fertile ground for growth of the MJI, such that the models evolve quickly into self-gravitational runaway, within 2 orbits. The evolution of model MS1 is similar to that of H1, in which perturbations grow primarily via swing amplification, although magnetic fields and stronger gravity expedite the instability in MS1.

In Figure 8$`a`$, we plot a snapshot of model MS3 at $`t/t_{\mathrm{orb}}=1.78`$ in the frame comoving with the spiral pattern (that is, the frame of the simulation). Surface density in logarithmic color scale shows the structure of spurs forming almost perpendicularly to the shock, with an average spacing of $`0.8\mathrm{kpc}`$. The velocity field clearly shows the background shear of galactic rotation across the box, the streaming motion of fluid elements within the arm, and the reversal of shear in the region of high density compression. Since it amounts to only $`5\%`$ of the amplitude of background velocity, however, the perturbed velocity associated with spur formation is difficult to see in this representation. The spurs move in the $`y`$-direction with $`v_y0.50R_0\mathrm{\Omega }_0`$ with respect to the global spiral pattern, which implies that they follow the background galactic rotation very closely (since $`v_{\mathrm{inertial}}=R_0\mathrm{\Omega }_p+v_y0.5R_0\mathrm{\Omega }_0+v_y`$). To see the shape of spurs clearly, we transform to a frame in which spurs remain stationary and the left $`x`$-boundary corresponds to the initial shock front location. The evolved density structure (color scale) and streamlines (dotted lines) in this frame are shown in Figure 8$`b`$. We also show the wavefront formed by spurs (solid line) drawn from equation (12), with $`k_x(0)/k_y=2.5`$ chosen to match the shape of spurs. Note that except in the region very close to the shock front, the streamlines in the stationary-spur frame follow the kinematic wavefronts of spurs fairly well.

Figure 9 displays snapshots of density and field configurations of models MS1 at $`t/t_{\mathrm{orb}}=4.0`$ and MS2 at $`t/t_{\mathrm{orb}}=2.1`$. The late-stage evolution of MS3 is similar to that of MS2. Our simulation resolution is insufficient to follow the evolution of these models further. When spurs reach high enough density, they begin to fragment along their length. Each of these fragments has an average mass of about 2.0% and 1.5% of the total for models MS1 and MS2, respectively, which correspond to $`4\times 10^6\mathrm{M}_{}`$. Notice that this clump mass is roughly equal to the local Jeans mass, $`M_{\mathrm{J},\mathrm{sp}}`$, at the spiral arm density peak

$$M_{\mathrm{J},\mathrm{sp}}\frac{c_s^4}{G^2\mathrm{\Sigma }_{\mathrm{max}}}=7\times 10^6\mathrm{M}_{}\frac{Q_{\mathrm{sp}}^2}{Q_0}\left(\frac{c_s}{7.0\mathrm{km}\mathrm{s}^1}\right)^3\left(\frac{\kappa _0}{36\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1}\right)^1.$$
(14)
From Table 1, one can see that $`M_{\mathrm{J},\mathrm{sp}}`$ ranges $`(2.54)\times 10^6\mathrm{M}_{}`$ for our magnetized models. Since the spur spacing is $`25\lambda _{\mathrm{J},\mathrm{sp}}`$, and the surface density has large gradients across the arm width of $`12\lambda _{\mathrm{J},\mathrm{sp}}`$, one should however not think of the clumps as being gathered isotropically from uniform density regions.

Some fragmentation occurs within the spiral arm (and fragments subsequently climb out of the potential well), while some fragments that form off the arm are rapidly carried away deep into the interarm region. Although the mass conversion efficiency in clustered star formation is uncertain, the collapsing fragments of our models may represent the entities that develop into prominent H II regions, as observed in both arm and interarm locations in spiral galaxies. Note that as Table 1 lists, the spacing along the arm of spurs in models MS2 and MS3 are about $`34`$ times the arm width, which is consistent with observed separation of H II regions and star complexes in spiral galaxies (Elmegreen & Elmegreen, 1983).

### 4.3 MHD Simulations with Equipartition Magnetic Fields

Figure 10 shows the initial density distributions and the evolution of maximum surface densities for MHD simulations with $`\beta _0=1`$ (i.e., $`v_\mathrm{A}=c_s`$). Since magnetic fields exert pressure as if they were adiabatic gas with an “effective” index of two for one-dimensional compression (cf. Shu 1992), they potentially can reduce the density compression greatly even in an isothermal shock. The maximum densities in the resulting steady-state shock profiles in models ME1$``$ME4 are less than 3, and the maximum to minimum density ratios are smaller than 6 as well. The minimum value of $`\beta `$ is $`0.30.5`$. With the corresponding values $`Q_{0,c}0.860.89`$ for the existence of equilibrium configurations, the chosen models having $`Q_0=1.21.6`$ are quite far from the margins of susceptibility for quasi-axisymmetric instability. Nevertheless, these $`\beta _0=1`$ models respond with strong amplification of perturbations that vary along the spiral arm, forming self-gravitating structures within $`2.5`$ orbits, as seen Figure 10$`b`$.

In order to study the morphology and kinematics of spurs, we display in Figure 11 various aspects of the surface density and velocity structure of model ME3 at $`t/t_{\mathrm{orb}}=1.27`$. Figure 11$`a`$ clearly shows that nine spurs, whose density is displayed in logarithmic color scale, grow out of the spiral arm. Velocity vectors measured in the spiral arm frame show the background shear, reversed shear, and streaming motions associated with passage through the spiral arm, but as for model MS3 (Fig. 8$`a`$) they do not readily evidence the perturbed velocity components that create spurs. The spurs of model ME3 have an average spacing of $`\lambda _y0.7\mathrm{kpc}`$, and move along the arm with $`v_y0.47R_0\mathrm{\Omega }_0`$. Figure 11$`b`$ shows that the shape of the spurs is in excellent agreement with the kinematic prediction for wavefronts (solid lines) based on background shear and expansion from equation (12) with $`k_x(0)/k_y=0.2.`$444With $`k_x(0)/k_y>0`$, the spurs are trailing when they leave the shock front ($`\tau =0`$). Since $`k_x(0)/k_y=0.2`$ is relatively small, however, any appreciable $`\tau `$ in eq. (12) gives $`k_x/k_y=𝒯<0`$ so that wavefronts become leading almost immediately downstream of the shock front, as Fig. 11$`b`$ shows. Convergence of streamlines (dotted lines) towards the regions of high density in spurs suggests that they grow by accumulating material mostly in the $`y`$-direction, which in turn implies that spurs in model ME3 are products of MJIs (Elmegreen, 1987a; Kim & Ostriker, 2001).

To provide a feeling for spur-region kinematics as it could be observed with a radio interferometer or Fabry-Perot instrument, we create sample synthetic maps. Imagine a spiral galaxy whose local disk can be identified with model ME3 at $`t/t_{\mathrm{orb}}=1.27`$. The galaxy is assumed to be inclined arbitrarily by 50 degrees with respect to the plane of sky. For fictitious observations, we select two viewing angles: parallel to the spiral arm and 40 degrees off the arm. In both cases, the target areas are rotating away from the observer. The resulting line-of-sight velocities (color scale) and “deprojected” column density (contours) are presented, respectively, in Figures 11$`c`$ and 11$`d`$. The numbers in the color bars are in the units of km s-1. When viewed along the spiral arm, the line-of-sight velocity is equal to $`v_y\mathrm{sin}\alpha `$ (where $`\alpha `$ is the angle of the disk plane with respect to the plane-of-sky). With this observer position angle, the signature of nonuniform streaming motions influenced by the external potential is evident, as manifested by the small change in color gradient across the arm in Figure 11$`c`$. The line-of-sight velocity has no discontinuities in this case. When the line-of-sight is oblique to the arm, on the other hand, the line-of-sight velocity is a mixture of the velocities both parallel and perpendicular to the arm. Correspondingly, the signature of parallel streaming motion in the line-of-sight velocity is suppressed relative to the signature of the large discontinuity in $`v_x`$ across the shock, as shown in Figure 11$`d`$ . A careful examination of Figures 11$`c`$ and 11$`d`$ reveals that the line-of-sight velocity is correlated with spurs, with high-density regions having negative gradients in the velocity ($`dv_y/dy<0`$), indicative of converging flow. The amplitude of the variation in the line-of-sight velocity along the arm is typically $`5`$ km s-1 in the spur regions, which is almost comparable to the variation of streaming motion (i.e., shear in $`x`$ of $`v_y`$) in the arm.

As mentioned before, the primary mechanism for structure formation in $`\beta _0=1`$ models is the MJI, in which threaded near-azimuthal magnetic fields transport angular momentum out of growing perturbations. Time evolution of the linear MJI in a disk with uniform density, uniform magnetic fields, and normal linear shear has been studied by Elmegreen (1987a), Gammie (1996), Fan & Lou (1997), and very recently by Kim & Ostriker (2001) – in which nonlinear evolution was also studied. In previous works, however, the effects of varying shear and background compression and expansion across a spiral shock were not taken into account. Because of the evidence from our simulations of the importance of MJIs, it is valuable to generalize the analytical work of B88 on hydrodynamical self-gravitating instabilities in spiral arms to include magnetic effects. This will enable us to make direct comparisons between the linear evolution of MJIs in non-uniform media and the results of our numerical simulations.

In Appendix A, we present the linearly perturbed MHD equations in the Lagrangian frame moving with the unperturbed background flow through a spiral arm. In order to integrate the governing set of equations, we need to prescribe an equilibrium density configuration. For an exemplary run, we adopt that of model ME3 (with $`Q_0=1.5`$, $`\beta _0=1`$, and $`F=3\%`$) as shown in Figure 3. Just as in B88, we follow the growth of a disturbance assumed to maintain phases of an initial locally plane-wave form $`e^{i(k_x(0)x_0+k_yy_0)}`$, where $`x_0`$ and $`y_0`$ are the initial position of a fluid element. We are free to choose $`k_x(0)/k_y`$ and $`K_y`$ of any input perturbations. Taking $`\delta \mathrm{\Sigma }/\mathrm{\Sigma }0`$ and $`\delta u=\delta v=dm=0`$ as initial conditions at the shock front, and varying $`k_x(0)/k_y`$ and $`K_y`$, we directly integrate equations (A7)$``$(A10) over time. For a given wavelet, we measure the total amplification factor and the time $`\tau _{\mathrm{grow}}`$ when maximum amplification occurs, and plot these values in Figure 12 as solid and dotted contours, respectively. Overall, modes with $`K_y<0.5`$ are favored for growth, in sharp contrast to MJIs in a uniform medium where dominant modes have $`K_y\frac{1}{2}\frac{3}{4}`$. This is because $`k_{\mathrm{J},\mathrm{sp}}`$ used in the definition of $`K_y`$ is in fact the maximum possible value in the medium. The value of $`k_\mathrm{J}=2\pi G\mathrm{\Sigma }/c_s^2`$ averaged over the arm is lower, rendering the “average” value of $`k_y/k_\mathrm{J}`$ larger than $`K_yk_y/k_{\mathrm{J},\mathrm{sp}}`$.

Since MJIs generally require the wavenumber of perturbations to be less than the local Jeans wavenumber for instantaneous growth, the total amplification factor depends on how long the kinematics of the background flow can keep $`k_x`$ small such that $`(k_{x}^{}{}_{}{}^{2}+k_{y}^{}{}_{}{}^{2})^{1/2}`$ is smaller than the local $`k_\mathrm{J}`$ (which itself decreases as the surface density drops outside the spiral arm). While $`k_x`$ is a linearly increasing function of time in a normal, uniform-shear case, the evolution of $`k_x`$ in the spiral arm is not generally monotonic, determined instead by equation (12). It turns out, as Figure 12 shows, that the MJIs associated with the spiral arm prefer wavelets that are initially slightly trailing at the shock front ($`k_x(0)/k_y>0`$). This is because the shear reversal inside the spiral arm produces an interval in which $`k_x`$ initially decreases, followed by the return of normal shear downstream which then increases $`k_x`$. The expansion effect also tends to reduce $`k_x`$ if it is initially positive. Starting a wavelet as trailing thereby extends the time span during which MJIs operate efficiently, and thus yields maximum amplification.

In Figure 12, we also indicate the parameters ($`K_y=0.45`$ and $`k_x(0)/k_y=0.2`$) determined from the numerical simulation for model ME3 (see Fig. 11) as a rectangular box. The ambiguities in determining $`k_y`$ from an integral number of spurs and in finding $`k_x(0)/k_y`$ that gives the best fit for the shape of spurs are represented by the box size. The corresponding amplification factor and growth time from the linear-theory integration are $`6080`$ and $`\tau _{\mathrm{grow}}3.13.5`$, respectively. While the largest predicted linear amplification available for the system is attained at the same $`k_x(0)/k_y=0.2`$, but at a larger length scale with $`K_y=0.32`$, this mode (having $`\tau _{\mathrm{grow}}=3.6`$) would take longer to achieve the maximum growth. It appears that the nonlinear system “compromises” between maximum amplification and earliest growth, in selecting the mode that dominates the evolved state of the simulation.

In Figure 13, we compare the variation in the perturbed surface density distribution within the spurs with the prediction of amplification from the linear analysis. In doing this, we use the equivalence between the spatial variable $`x`$ in the simulation and the Lagrangian distance traveled downstream by a fluid element in the background state, in the linear analysis. Linear evolution of the perturbed surface density with $`K_y=0.45`$ and $`k_x(0)/k_y=0.2`$ is plotted as a function of the dimensionless Lagrangian time variable $`\tau `$, and the corresponding downstream distance from the shock front. Note that $`\tau `$ varies from 0 to $`2\pi m^1(1\mathrm{\Omega }_p/\mathrm{\Omega }_0)^1`$ for $`x/L_x=01`$, corresponding to $`\tau =02\pi `$ for two-armed spirals with $`\mathrm{\Omega }_p=\mathrm{\Omega }_0/2`$. We select three prominent spurs in Figure 11$`b`$ and plot with various curves their perturbed surface densities measured along their length. Although densities near the shock front fluctuate with small amplitudes due to the interactions with spurs’ trailing tails, the overall density distribution and the positions where spurs achieve maximum densities are in remarkably good agreement with the results of linear analysis.

The growth of perturbations depends sensitively on $`Q_0`$. For instance, model ME4 with $`Q_0=1.6`$ takes about one more orbital time than model ME3 with $`Q_0=1.5`$ to be fully nonlinear, although initial density distributions are quite similar to each other (see Fig. 10). This is not because the linear analysis for $`Q_0=1.6`$ predicts a larger value of $`\tau _{\mathrm{grow}}`$ than for $`Q_0=1.5`$, but because the perturbations keep growing as they pass successively through spiral arms. As $`Q_0`$ increases further, linear amplification factors become smaller, so that higher-$`Q_0`$ models need even longer time for development. We found that a model with $`Q_0=1.7,\beta _0=1`$, and $`F=3\%`$ (not listed in Table 1) becomes only moderately nonlinear at $`t/t_{\mathrm{orb}}=4`$, while models with $`Q_01.8`$ and the same $`\beta _0`$ and $`F`$ remain linear until the end of simulations ($`t/t_{\mathrm{orb}}=5`$).

Finally, Figure 14 displays developing and fragmenting spurs for models ME3 at $`t/t_{\mathrm{orb}}=1.59`$ and ME4 at $`t/t_{\mathrm{orb}}=2.55`$. The separation between two neighboring spurs is about $`0.70.8`$ kpc, which is very close to the most unstable Jeans wavelength at the spiral density peak ($`K_{y,\mathrm{max}}0.45`$; see Table 1), showing that the stabilizing contribution from epicyclic motions is almost negligible when the magnetic field is strong. Compared to $`\beta _0=10`$ models, spurs in the $`\beta _0=1`$ models are stronger and extend farther away from the spiral arm. Fragmentation of spurs when the magnetic field is stronger thus occurs mostly in the interarm region. For these runs, we find that each collapsing fragment has an average mass of about 1.5% ($`4\times 10^6\mathrm{M}_{}`$) of the total mass, which is again very close to the local Jeans mass $`M_{\mathrm{J},\mathrm{sp}}`$ at the peak density of spiral arms (see eq.  and Table 1). One cannot resist speculating that bright, interarm H II regions, as seen for example in the optical image of M51, could have originated as condensations in spur structures very like those found by our simulations and analyses.

## 5 Discussion

### 5.1 Summary of Modeling and Results

In this paper, we have investigated the dynamical interaction of a magnetized, 2D shearing flow – representing a local patch of the height-integrated ISM – with an external gravitational potential – representing a local portion of a spiral stellar arm. Our primary interest was in exploring the ways in which self-gravity leads to growth of intermediate-scale structure, and how these intermediate-scale structures subsequently fragment. We show that in magnetized systems, the characteristic structures that grow are spur-like features that jut out at regular intervals from the spiral arms; within these spurs, dense condensations may grow and travel out into the interarm region. The spur structures in our models bear remarkable resemblance to the conspicuous features branching out of spiral arms in recent high-resolution HST images of the Whirlpool galaxy, and also apparent (less spectacularly) in other grand-design spirals.

The technical centerpiece of our work is a set of time-dependent, numerical MHD simulations. The model disks are infinitesimally thin and maintain a constant isothermal sound speed $`c_s`$. In the absence of stellar spiral-arm forcing, the disks are assumed to have uniform velocity shear, uniform surface density characterized by the Toomre $`Q_0`$ parameter (see eq. ), and uniform near-azimuthal magnetic fields, with $`\beta _0c_s^2/v_\mathrm{A}^2`$ measuring the mean field strength at the disk midplane (see eq. ). We do not allow for spacetime variations of the magnetic scale height. In addition to its own self-gravity, the gaseous flow is subject to an explicitly-introduced sinusoidal external potential (see eqs.  and ), modeling a local stellar spiral arm. The back reaction of stars to the gravitational field of gas is not, however, included in the present study. The spiral arm is assumed to be tightly wound and rigidly rotating at half the local orbital rate. We evolve the ideal MHD equations in the local frame comoving with the stellar arm, with two orthogonal axes corresponding to the directions perpendicular and parallel to the arm (see Roberts 1969).

Imposing translational symmetry along the arm, we first construct one-dimensional equilibrium spiral shock configurations and examine their stability to axisymmetric perturbations. Even if the spiral perturbation is only a tiny fraction of the background axisymmetric field, induced radial velocities are generally supersonic, easily forming a shock front. We find that larger amplitude background spiral arms (larger $`F`$) result in stronger shocks, while the shock-forming tendency of gas decreases as the relative importance of self-gravity to spiral forcing increases (cf. Lubow, Balbus, & Cowie 1986). Self-gravity enhances the arm-interarm density contrast, while magnetic pressure reduces it. When $`Q_0`$ is sufficiently small, no equilibrium profiles can be found from our time-dependent integrations, because either the structures become gravitationally unstable (when magnetic fields and/or stellar spiral perturbations are weak), or because the given background conditions and adopted spiral pattern speed do not allow stationary shocks to exist (when fields and/or spiral arms are strong). Critical values of $`Q`$ at the density peak are found to be $`Q_{\mathrm{sp}}0.8,`$ 0.5, and 0.4 for $`\beta _0=\mathrm{}`$, 10, and 1, respectively.

For our two-dimensional simulations, we start with one-dimensional equilibrium density profiles calculated as above, apply low-amplitude white-noise perturbations, and monitor their nonlinear growth. We find that gaseous spurs naturally form as a consequence of gravitational instabilities inside spiral arms with input modes having $`𝐤`$ nearly along the spiral arms growing preferentially. The chief physical mechanism for spur formation is the magneto-Jeans instability, in which the magnetic tension force breaks the loop of stabilizing epicyclic motions (Elmegreen, 1987a; Kim & Ostriker, 2001). Although swing amplification can also play a role in forming spur structures when the magnetic field is weak and the arm-interarm density contrast is moderate, purely hydrodynamic systems that are stable to quasi-axisymmetric perturbations are generally also found to be stable to nonaxisymmetric perturbations. We thus suggest that magnetic effects are key for producing the gaseous spurs observed in real galaxies.

When the mean magnetic field is in thermal equipartition ($`\beta _0=1`$), we find that the separation of spurs is more-or-less consistent with the most unstable Jeans wavenumber at the density peak, while with a sub-equipartition magnetic field ($`\beta _0=10`$), the formation of spurs requires larger-scale perturbations to enhance self-gravity (see Table 1). The characteristic scale of spurs formed in our simulations is about $`\lambda _y2.5\lambda _{\mathrm{J},\mathrm{sp}}`$ (or $`750`$ pc for our model parameters) for the models that are unstable to the magneto-Jeans instability, where $`\lambda _{\mathrm{J},\mathrm{sp}}`$ is the local Jeans wavelength at the spiral-arm density peak. The corresponding ratio of the spacing to the arm thickness is $`\lambda _y/W3.5`$ and 1.5, respectively, for $`\beta _0=10`$ and 1 models. These ratios are consistent with the observational results for the distribution of H II regions along spiral arms (Elmegreen & Elmegreen, 1983).

The shape of the spurs sculpted by the shearing and expanding background flow in our models is in excellent agreement with equation (12), originally derived by B88. Reversed shear within the strongly compressed region tends to decrease the $`x`$-wavenumber of the perturbations, so that slightly trailing waves entering the spiral shock are preferred as having the longest growing phase under MJIs. Entering waves first shear around and expand to create nearly-radial wavefronts as they exit the arm. Far into the interarm region, where shear returns to the “normal” (slower-outside) sense, the wavefronts become strongly trailing.

In the nonlinear stage of evolution, spurs experience fragmentation to form gravitationally bound clumps, the typical mass of which corresponds roughly to the local Jeans mass at the density peak inside spiral arms (or $`4\times 10^6\mathrm{M}_{}`$ for our model parameters). The positions where the fragmentation occurs include interarm regions as well as within spiral arms themselves. We propose that interarm condensations formed in this way could evolve into bright H II regions in interarm regions.

### 5.2 Application to Spiral Galaxies

Now we apply our simulation results to a grand-design spiral, M51. We are particularly interested in a section of the southern spiral arm at 2.7 kpc from the center, where spurs can be relatively easily identified and observed data are available. Here, we adopted a galactic distance of 9.6 Mpc to M51 (Sandage & Tammann, 1975). From the recent combined WFPC2/NOAO Hubble Heritage image of M51 (Scoville et al, 2001), we recognize five spurs in this portion of the arm with an average spacing of $`\lambda _y350`$ pc. From radio interferometry and single dish observations, $`\mathrm{\Omega }_078\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1`$, $`\mathrm{\Sigma }_0=80\mathrm{M}_{}\mathrm{pc}^2`$, and $`\mathrm{\Sigma }_{\mathrm{max}}=170\mathrm{M}_{}\mathrm{pc}^2`$ (Rand, 1993; Garcá-Burillo et al, 1993). With $`c_s10\mathrm{km}\mathrm{s}^1`$, the values listed above give $`\lambda _{\mathrm{J},\mathrm{sp}}140`$ pc, corresponding to $`K_{y,\mathrm{max}}\lambda _{\mathrm{J},\mathrm{sp}}/\lambda _y0.4`$. This is consistent with Table 1, which suggests that spurs produced by MJIs have $`K_{y,\mathrm{max}}0.30.5`$.

Giant H II regions are well distributed along the dominant spiral arms in grand design spirals, but a significant fraction of them are also found in deep interarm regions. Scoville et al (2001) reported that about 55% of identified H II regions are in interarm regions, although overlapping of H II regions in the HST image of M51 could reduce their numbers in the arms. Given the nominally low surface density and high shear conditions in interarm regions, the presence of numerous interarm H II regions requires an explanation. From our model simulations, we suggest that massive self-gravitating clumps resulting from nonlinear fragmentation of spurs can develop in both arm and interarm regions. This implies that interarm star formation may be engendered by the same dynamical process, namely large-scale gravitational instabilities initiated in spiral arms, as arm star formation. Indeed, there is no significant observational difference in the sizes and electron densities of the H II regions between spiral arm and interarm regions (Scoville et al, 2001).

The fact that the shape of spurs, reflecting the background flow characteristics, is well described by equation (12) may provide an independent way to determine global spiral pattern speeds. The pitch angle, angular velocity of galactic rotation, and epicyclic frequency are observable quantities, and $`k_y`$ can be determined from the mean separation of observed spurs. Provided that the surface density distribution is known, therefore, $`k_x(0)`$ and $`\mathrm{\Omega }_p`$ (see eq. ) can be determined simultaneously by finding a best fit of $`𝒯`$ to the observed overall shape of spurs. If we can further constrain $`k_x(0)/k_y`$ for example from the linear theory as the dominant mode, or from high resolution observations that allow determination of the slope of spurs near the density maximum ($`\tau 0`$), then $`\mathrm{\Omega }_p`$ can be uniquely determined. Conversely, if we know the pattern speed by other means, equations (12) and (13) can be used to predict the variation of surface density in spiral arms, as a consistency check of direct determinations.

For model simulations presented in this paper, we took a patch of a galactic disk inside the corotation radius of a steady spiral pattern, and showed that spurs jut outwards from the arms and become increasingly trailing at larger radii. We have also run spiral-arm interaction simulations for regions outside the pattern’s corotation (not listed in Table 1) and found that the induced spurs jut inwards from the local arm segments and become trailing at smaller radii. These morphological differences in spurs between inside and outside corotation could help to locate corotation radii of spiral patterns in spiral galaxies. The present lack of evidence for spurs on the insides of arms (B. Elmegreen, personal communication) could be because the arms are generally relatively weak outside corotation, or perhaps because resolution is insufficient. Of course, it is possible that when other galaxies are observed at resolution comparable to the recent M51 image, and when larger areas are surveyed, inward-jutting spurs might be detected.

### 5.3 Outstanding Issues

In this paper, the 2D simulations of purely hydrodynamic flow passing through spiral arms show comparatively stable behavior, largely because the hydrodynamic models chosen for these runs have relatively large $`Q_{\mathrm{sp}}>0.9`$ at the density peak. Balbus & Cowie (1985) argued that for flows with pressure-density relation $`\mathrm{\Pi }\mathrm{\Sigma }^\gamma `$, the values of the Toomre parameter at the point of maximum compression, $`Q_{\mathrm{sp}}Q_0\gamma ^{1/2}(\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_0)^{\gamma /21}`$, is the most important parameter for characterizing the system’s gravitational response (here $`\mathrm{\Sigma }_0`$ and $`Q_0`$ are the equivalent parameters for a uniform disk). Since $`Q_{\mathrm{sp}}`$ is the smallest local value in a non-uniform medium, $`Q_{\mathrm{sp}}<1`$ does not always ensure quasi-axisymmetric instability will occur, because for local instability $`Q`$ must be smaller than unity over at least a radial distance $`\lambda _{\mathrm{max}}=2c_s^2/G\mathrm{\Sigma }=2\gamma c_{s,0}^2(G\mathrm{\Sigma }_0)^1(\mathrm{\Sigma }/\mathrm{\Sigma }_0)^{\gamma 2}`$. Our results for quasi-axisymmetric stability indicate that critical $`Q_{\mathrm{sp}}`$ values are $`0.730.98`$ for unmagnetized cases (see §3 and Fig. 4). It seems that at least one additional parameter, $`W`$, which describes arm width (or the gradient in the background density profile), combines with $`Q_{\mathrm{sp}}`$ to determine the system’s response to two-dimensional (spur-forming) perturbations. For example, models H2 and H3 having smaller $`Q_{\mathrm{sp}}`$ remain stable (with fluctuating density), while model H1 with larger $`Q_{\mathrm{sp}}`$ (but smaller $`Q_0`$) is subject to (mild) swing amplification. For a given $`Q_{\mathrm{sp}}`$, models with larger $`W`$ – and thus a smaller arm-to-interarm density contrast and $`Q_0`$ – are usually more susceptible to swing amplification. In part, this is because excessive compression in the spiral arm reverses the normal direction of shear (see eq. ), which tends to suppress the swing mechanism. For wavelets that enter the spiral shock with leading orientation ($`k_x<0`$), B88 showed that the tendency for “normal” rotation of leading wavefronts driven by expansion (see eq. \[3.15\] of B88) can compensate somewhat for locally-reversed shear. We find that for cases where quasi-axisymmetric modes are stable, however, amplification from “natural” perturbations (mostly trailing) does not produce strong spurs in unmagnetized systems.

From B88’s results with $`\gamma =0.5`$, $`Q_{\mathrm{sp}}<0.7`$ appears to be required in order for nonaxisymmetric perturbations to grow significantly enough to form nonlinear structures. With a small value of $`\gamma `$, moderate compression factors easily produce $`Q_{\mathrm{sp}}`$ as low as 0.6 from mean ISM conditions. However, quasi-axisymmetric modes may be unstable as well when $`Q_{\mathrm{sp}}`$ is very small. By adopting an isothermal equation of state, we obtain larger values of $`Q_{\mathrm{sp}}`$ for a given compression factor than with a softer pressure-density relation. In addition, because our 2D models are initiated with 1D profiles that are quasi-axisymmetrically stable, we were limited to a range of $`Q_{\mathrm{sp}}>0.9`$ for our hydrodynamic models. Although introduced to some extent to achieve numerical stability in our computations, the requirement of quasi-axisymmetric stability is likely to appear in nature as well. Since we observe gaseous spiral arms and associated dust lanes in spiral galaxies, the compressed gaseous regions are probably quasi-axisymmetrically stable; otherwise, these arms would not be long-lived – compressed gas would turn quickly into stars.

The question, then, of whether spurs could form from purely hydrodynamic (swing-like) instabilities within arms is reduced to finding a self-consistent density distribution that is quasi-axisymmetrically stable yet with $`Q_{\mathrm{sp}}`$ small enough for 2D instabilities to grow. The real ISM is multiphase and turbulent, in which hot, warm, and cold phases coexist (e.g., Field, Goldsmith, & Habing 1969; Cox & Smith 1974; McKee & Ostriker 1977; Heiles 2001). Although optically thin gas can be described by polytropic relations between density and thermal pressure, this is only in a piecewise sense with a fairly uncertain $`\gamma `$ (see, e.g., Vázquez-Semadeni et al 2000). The turbulent motions of cold clouds, although often treated as an effective pressure, may actually have a rather different dynamical response to compression. A realistic assessment of hydrodynamic instabilities on spur and cloud formation, therefore, will require a more accurate treatment of the small-scale thermal and dynamical properties of the medium. An important direction for future research, thus, will focus on determining whether better treatment of “microphysics” renders the purely-hydrodynamic (swing) mechanism competitive with the MHD mechanism (MJI) for forming spurs and condensing clouds.

In the present paper, we have studied the formation and fragmentation of spurs due to self-gravitating instabilities within spiral arms, with magnetic fields parallel to the galactic plane abetting the process. Since they are restricted to a thin-disk geometry, the present models do not capture the potential dynamical consequences of magnetorotational instabilities (MRIs) and the Parker instability in themselves or in connection to MJIs. MRIs (e.g., Balbus & Hawley 1998), which exist only in 3D systems, are known to generate MHD turbulence, which may make an important contribution to the amplitude of ISM random motions in galactic disks (Sellwood and Balbus, 1999) as well as exciting density perturbations, and may also provide angular momentum transport within disks – possibly affecting the growth of condensations.

The Parker (1966) instability has long been thought to be important in the formation of OB associations, giant H II regions, and giant molecular clouds along spiral arms (e.g., Mouschovias, Shu, & Woodward 1974; Blitz & Shu 1980), because it is able to grow condensations at lateral wavelengths comparable to the disk scale height – as required for clouds of mass $`10^5\mathrm{M}_{}`$. Recent numerical studies indicate that the Parker instability alone cannot be the main formation mechanism for giant clouds in general galactic disks, because the column density enhancement is only a factor of three (Kim et al 2000; Santillán et al 2000; see also Elmegreen 1995). However, the coupling of the Parker instability with other dynamical processes could produce significant effects on structure formation in galaxies. Basu, Mouschovias, & Paleologou (1997) suggested that nonlinear triggering of the Parker instability by spiral shocks could enhance density enough to form giant cloud complexes, and the nonlinear simulations of Chou et al (2000) (neglecting shear) evidence this sort of behavior. Extension of the current work into three dimensions has the potential to show how “seeds” sown by the Parker and/or MRI modes in spiral arms may grow and condense into self-gravitating giant molecular clouds.

It is a pleasure to acknowledge stimulating conversations and communications with B. Elmegreen, J. Stone, S. Vogel, F. Shu, N. Scoville, and K. Sheth. We are also grateful to an anonymous referee for constructive comments. This work was supported by NASA grants NAG 53840 and NAG 59167.

## Appendix A Linear Analysis

In this appendix, we provide the perturbed equations to be integrated to obtain the linear-theory solutions shown in Figures 12 and 13. These solutions represent examples of the evolution of low-amplitude, nonaxisymmetric perturbations (specifically, initial local plane waves) in a magnetized, non-uniform disk, subject to a background flow containing compression and expansion, and varying shear, as imposed by the self-consistent equilibrium response to an external spiral potential. Local dynamical instabilities in unmagnetized disks with similar background flow properties have been studied by Balbus & Cowie (1985) for quasi-axisymmetric modes and by B88 for two-dimensional perturbations. Here, we extend B88’s work to include the effect of magnetic fields.

We begin by considering a steady-state equilibrium configuration represented by density $`\mathrm{\Sigma }`$, velocity $`𝐯=u\widehat{𝐱}+v\widehat{𝐲}`$, and magnetic field $`𝐁=B_y\widehat{𝐲}`$ (see for example §3), in the spiral-arm coordinates built in §2.1. The equilibrium is quasi-axisymmetric in the sense that flow quantities vary only along $`\widehat{𝐱}`$, the coordinate axis perpendicular to the arm, which lies at an angle $`i`$ with respect to the local radial direction. Imposing small amplitude Eulerian perturbations represented by $`\delta `$, we linearize equations (1)-(5):

$$\frac{d}{dt}\left(\frac{\delta \mathrm{\Sigma }}{\mathrm{\Sigma }}\right)=\frac{\delta u}{x}\frac{\delta v}{y}\delta u\frac{d\mathrm{ln}\mathrm{\Sigma }}{dx},$$
(A1)
$$\frac{d\delta u}{dt}=\delta u\frac{du}{dx}+2\mathrm{\Omega }_0\delta v\frac{}{x}\left(c_s^2\frac{\delta \mathrm{\Sigma }}{\mathrm{\Sigma }}+\delta \mathrm{\Phi }_g\right)+\frac{B_y}{4\pi \mathrm{\Sigma }}^2\delta m+\frac{1}{4\pi \mathrm{\Sigma }}\frac{dB_y}{dx}\left(\frac{\delta m}{x}+B_y\frac{\delta \mathrm{\Sigma }}{\mathrm{\Sigma }}\right),$$
(A2)
$$\frac{d\delta v}{dt}=\left(\frac{\kappa _0^2}{2\mathrm{\Omega }_0}+\frac{dv}{dx}\right)\delta u\frac{}{y}\left(c_s^2\frac{\delta \mathrm{\Sigma }}{\mathrm{\Sigma }}+\delta \mathrm{\Phi }_g\right)+\frac{1}{4\pi \mathrm{\Sigma }}\frac{dB_y}{dx}\frac{\delta m}{y},$$
(A3)
$$\frac{d\delta m}{dt}=B_y\delta u,$$
(A4)
$$^2\delta \mathrm{\Phi }_g=4\pi G\delta (z)\delta \mathrm{\Sigma },$$
(A5)
where the Lagrangian time derivative following the background flow is denoted by

$$\frac{d}{dt}\frac{}{t}+u_T\frac{}{x}+v_T\frac{}{y},$$
(A6)
and the perturbed vector potential $`\delta m`$ is defined through $`\delta 𝐁\times (\delta m\widehat{𝐳})`$.

Following B88, we now consider a Lagrangian frame comoving with the background flow $`𝐯_T`$, allowing for the local normal or reversed shear and expansion velocity fields. The frame is initially located at the postshock density peak. We adopt the ansatz that an applied plane-wave disturbance locally preserves sinusoidal variations, so that the stability of a local patch can be assessed by following the temporal evolution of wavelets of arbitrary initial local $`𝐤(0)`$. Strictly speaking, this sort of analysis should be applied only when the initial wavelength is small compared to the width of the spiral arm. Formally, however, we can choose $`𝐤(0)`$ freely. Since the background variables are independent of $`y`$, the $`y`$-wavenumber $`k_y`$ of the perturbations remains fixed throughout the linear evolution, but the $`x`$-wavenumber $`k_x`$ changes in response to the background flow. In particular, normal/reversed shear tends to increase/decrease $`k_x`$, and expansion/contraction tends to decrease/increase $`|k_x|`$. The net result of these various effects on the local value of $`k_x`$ was derived by B88; our equation (12) reproduces his formula for the local value of $`k_x/k_y=𝒯`$ as a function of the elapsed time since entering the shock (or, correspondingly the perpendicular $`x`$-distance traversed).

For consistency with other “local” simplifications we have made, we shall assume that spatial variations of the perturbed variables are much more rapid than those of the background state, so that we can ignore the last term in each of equations (A1)$``$(A3). The validity of this approximation can easily be demonstrated by computing $`rk_y(d\mathrm{ln}\mathrm{\Sigma }/dx)^1`$; if $`r1`$, we can safely ignore the terms containing the gradients of $`\mathrm{\Sigma }`$ and $`B_y`$ in equations (A1)$``$(A3) (recall that $`\mathrm{\Sigma }B_y`$ from flux conservation). Since $`d\mathrm{ln}\mathrm{\Sigma }/dxW^1`$, where $`W`$ is the arm width (cf. Table 1), it follows that

$$rWk_y33\left(\frac{K_y}{Q_0}\right)\left(\frac{\mathrm{\Sigma }_{\mathrm{max}}}{\mathrm{\Sigma }_0}\right)\left(\frac{W}{L_x}\right),$$

where $`K_yk_y/k_{\mathrm{J},\mathrm{sp}}`$, with the Jeans wavenumber $`k_{\mathrm{J},\mathrm{sp}}`$ at the density peak defined by equation (11). Using the background parameters listed in Table 1 with the identification of $`K_y=K_{y,\mathrm{max}}`$, corresponding to the most unstable modes, we find $`r45`$ for the $`\beta _0=1`$ models and $`r12`$ for the $`\beta _0=10`$ models. This implies that the WKB approximation we make here is somewhat marginal for the $`\beta _0=1`$ models, and probably breaks down for the $`\beta _0=10`$ models. For this reason, we shall apply the results of the linear analysis only for the $`\beta _0=1`$ models, particularly model ME3.

Defining the dimensionless variables, $`\tau t\mathrm{\Omega }_0`$, $`\delta \sigma \delta \mathrm{\Sigma }/\mathrm{\Sigma }`$, $`\delta \stackrel{~}{u}i\delta uk_{\mathrm{J},\mathrm{sp}}/\mathrm{\Omega }_0`$, $`\delta \stackrel{~}{v}i\delta vk_{\mathrm{J},\mathrm{sp}}/\mathrm{\Omega }_0`$, and $`\delta \stackrel{~}{m}i\delta mk_{\mathrm{J},\mathrm{sp}}/(B_y)`$, and applying the WKB approximation, we now rewrite equations (A1)-(A5) in dimensionless form (omitting the tilde) as

$$\frac{d\delta \sigma }{d\tau }=K_y(𝒯\delta u\delta v),$$
(A7)
$$\frac{1}{}\frac{d(\delta u)}{d\tau }=2\delta v\alpha K_y𝒯\left[1\frac{1}{|K_y|(1+𝒯^2)^{1/2}}\right]\delta \sigma \alpha \beta _0^1\sigma _{\mathrm{max}}K_y^2(1+𝒯^2)\delta m,$$
(A8)
$$\frac{d\delta v}{d\tau }=(2q_0)\sigma _{\mathrm{max}}\frac{\delta u}{}+\alpha K_y\left[1\frac{1}{|K_y|(1+𝒯^2)^{1/2}}\right]\delta \sigma ,$$
(A9)
$$\frac{d\delta m}{d\tau }=\frac{\delta u}{},$$
(A10)
where $`\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }`$, $`\sigma _{\mathrm{max}}\mathrm{\Sigma }_{\mathrm{max}}/\mathrm{\Sigma }_0`$, and $`\alpha (c_sk_{\mathrm{J},\mathrm{sp}}/\mathrm{\Omega }_0)^2=8(2q_0)\sigma _{\mathrm{max}}^2Q_0^2`$.

For temporal integrations of equations (A7)$``$(A10), we need to specify $`=(\tau )`$ and $`\sigma _{\mathrm{max}}`$ from a given background flow model: we adopt those of model ME3 shown in Figure 3. For initial conditions, we take $`\delta \sigma =1`$, and $`\delta u=\delta v=\delta m=0`$. We vary $`K_y`$ and $`k_x(0)/k_y`$ and compute amplification factors. The results are presented in Figures 12 and 13.