WIGGLE INSTABILITY OF GALACTIC SPIRAL SHOCKS: EFFECTS OF MAGNETIC FIELDS

Stellar spiral arms in disk galaxies greatly affect galaxy evolution in various ways (e.g., Buta & Combes 1996; Kormendy & Kennicutt 2004; Buta 2013; Sellwood 2014). They not only trigger and/or organize star formation by compressing gas into spiral shocks (e.g., Roberts 1969; Roberts & Yuan 1970; Shu et al. 1972, 1973) but also drive secular galaxy evolution (e.g., Lin & Shu 1964, 1966; Toomre 1964; Elmegreen 1995; Bertin & Lin 1996; Foyle et al. 2010). Observations commonly indicate that grand-design spiral arms abound with secondary structures called stellar "spurs" or gaseous/dust "feathers" that stretch out almost perpendicularly from the main arms into the interarm regions in a trailing configuration (e.g., Baade 1963; Lynds 1970; Weaver 1970; Elmegreen & Elmegreen 1983; Scoville & Rector 2001; Scoville et al. 2001; Kennicutt 2004; Willner et al. 2004; La Vigne et al. 2006; Corder et al. 2008; Silva-Villa & Larsen 2012; Schinnerer et al. 2013). While spurs and feathers have similar pitch angles and thus are thought to share the common origin (Elmegreen 1980; see also Puerari et al. 2014), what actually forms them has been a matter of considerable debate.

Existing theories for the formation of secondary structures differ in the relative role played by gaseous self-gravity to other hydrodynamic processes. Balbus & Cowie (1985) and Balbus (1988) studied a linear stability of spiral shocks to axisymmetric and non-axisymmetric modes, respectively, and showed that these arm substructures may form as a consequence of swing amplification of local density perturbations. Elmegreen (1994) showed that local perturbations can grow faster in the presence of magnetic fields that can remove the constraint of angular momentum conservation. These linear-theory predictions were confirmed by numerical simulations (Kim & Ostriker 2002, 2006; Shetty & Ostriker 2006) that showed that both magnetic fields and self-gravity are essential to form secondary structures resembling feathers in the nonlinear stage. Recently, Lee & Shu (2012) treated a stability of self-gravitating, magnetized spiral shock as a global, rather than local, problem in the direction perpendicular to the arm. They neglected galactic shear for perturbations while keeping non-inertial terms, and found semi-analytically that such shocks are prone to feather formation. Although these feather-forming instabilities are termed differently as the azimuthal instability, magneto-Jeans instability, and feathering instability by Elmegreen (1994), Kim & Ostriker (2002), and Lee & Shu (2012), respectively, they refer to the same gravitational instability of a rotating medium in which embedded magnetic fields play a destabilizing role. More recently, Lee (2014) extended the work of Lee & Shu (2012) to explore the parameter space of the feathering instability by varying the strength of magnetic fields and self-gravity.

In contrast, other studies argued that gaseous self-gravity may not be indispensable for the formation of feathers (e.g., Johns & Nelson 1986; Wada & Koda 2004; Dobbs & Bonnell 2006, 2007). These authors showed that secondary structures develop even in the absence of self-gravity, with their spacings smaller than those resulting from the self-gravitating instabilities mentioned above. In particular, Wada & Koda (2004) observed that small clumps form in strongly shocked layers that may grow into feathers, and termed the clump-forming mechanism the wiggle instability (WI). The WI appears prevailing in recent numerical simulations of spiral galaxies (Kim & Kim 2014) and even barred galaxies with strong dust-lane shocks (Kim et al. 2012a, 2012b; Kim & Stone 2012; Seo & Kim 2013, 2014). Based on the Richardson number criterion (e.g., Chandrasekhar 1961), they further proposed that the WI is originated from the Kelvin–Helmholtz instability (KHI) of a shear layer behind the shock. However, the linear analysis of Dwarkadas & Balbus (1996) showed that postshock flows in the presence of rapidly varying shear are stable to the KHI. Interestingly, Hanawa & Kikuchi (2012) argued that the WI may result spuriously from the difficulty in resolving a shock inclined to numerical grids.

To address the nature of the WI, Kim et al. (2014, hereafter Paper I) adopted the method of Lee & Shu (2012) and performed a linear stability analysis of non-self-gravitating, unmagnetized spiral shocks. Paper I found that the WI is physical rather than numerical, arising from the generation of potential vorticity (PV) at a distorted shock front, known as Crocco's theorem. Since gas in galaxy rotation passes through spiral shocks twice in every rotation for two-armed spirals, the PV contained in entropy-vortex waves can increase successively, in a Lagrangian sense, in the course of galaxy rotation. This sets up a Eulerian overstable normal mode along the azimuthal direction that grows exponentially. This is quite distinct from the KHI that relies on shear in the background flows. By relying on the growth of incompressible entropy-vortex modes, the WI is also different from the magneto-Jeans or feathering instability mentioned above that utilizes compressive acoustic modes. Paper I also confirmed that the growth rate and wavelength of the most unstable mode found by the linear stability analysis are consistent with the results of direct numerical simulations.

While the results of Paper I are informative to understand the nature of the WI, they are based on the models that do not consider magnetic fields pervasive in the interstellar medium (ISM). It has been well recognized that disk galaxies have a large-scale, ordered component as well as a small-scale, random component of magnetic fields (e.g., Wielebinski & Krause 1993; Heiles 1995; Beck et al. 1996; Beck 2001). In most disk galaxies, magnetic field directions based on polarized synchrotron radiation follow optical spiral structures fairly well, with their strength and pitch angle ranging typically $\sim 4-20\;\mu {\rm{G}}$ and $\sim 8^\circ -37^\circ$, respectively (see, e.g., Neininger 1992; Beck et al. 1996; Van Eck et al. 2015). Although there is an indication that dust lanes upstream of optical spiral arms have the strongest magnetic fields, some regular fields also extend well into interarm regions. It appears that small scale activities such as star formation and supernova explosions inside spiral arms disrupt regular magnetic fields and turn them into turbulent ones. When compressed by spiral shocks, these turbulent fields become as strong as, and sometimes even stronger than, the regular component inside the arms (e.g., Fletcher et al. 2011; Houde et al. 2013; Shneider et al. 2014). For gas surface densities of ${\rm{\Sigma }}\sim 10-100{M}_{\odot }\;{\mathrm{pc}}^{-2}$ inside the arms, these fields are close to equipartition strength with the thermal and turbulence energies, indicating that they are important in the gas dynamics associated with the arms (e.g., Beck et al. 1996; Chyży et al. 2000; Beck 2007).

In this paper, we extend Paper I by including the effects of magnetic fields on the stability of spiral shocks. As in Paper I, we consider a local shearing sheet of an infinitesimally thin gaseous disk that is assumed to be isothermal and non-self-gravitating. To simplify the situation further, we consider only the regular component of magnetic fields that is initially parallel to the imposed stellar spiral arms. One-dimensional (1D) solutions of equilibrium magnetized spiral shocks were already obtained by several authors (e.g., Roberts & Yuan 1970; Kim & Ostriker 2002; Lee 2014). By imposing two-dimensional (2D) perturbations onto such shocks, we will show that the presence of magnetic fields substantially suppresses the development of the WI mainly by reducing the density compression factor of a background spiral shock. The tension and pressure forces from the perturbed magnetic fields will turn out to have a minor effect on the WI. The stabilizing effect of magnetic fields is qualitatively consistent with the results of the previous MHD simulations (Kim & Ostriker 2006; Shetty & Ostriker 2006; Dobbs & Price 2008). We also run local MHD simulations to confirm the results of our linear stability analysis, and compare the pitch angles of the nonlinear structures formed in the simulations with the values reported by Puerari et al. (2014) for observed feathers in the grand-design spiral galaxies M51 and M81.

This paper is organized as follows. In Section 2, we introduce the basic MHD equations and our choices of the parameters, as well as the equilibrium shock solutions we use as a background state of the WI. In Section 3, we derive the linearized perturbation equations for our normal-mode stability analysis and present the shock jump conditions. We also give the computational method to obtain the dispersion relations. In Section 4, we present the resulting dispersion relations for 1D and 2D modes as well as the results of the MHD simulations, and discuss the effect of magnetic fields on the WI. In Section 5, we summarize our findings and discuss their astronomical implications.

2.1. Basic Equations

We consider magnetized spiral shocks in a galactic disk and study their stability to small-amplitude perturbations. The disk is assumed to be infinitesimally thin, isothermal with the speed of sound ${c}_{{\rm{s}}}$, and non-self-gravitating: the effect of self-gravity will be studied in a separate work.

Following Roberts (1969), we employ a local Cartesian coordinate system $(x,y)$ with the x- and y-axes representing the directions perpendicular and parallel to a local segment of a spiral arm, respectively (see also Shu et al. 1973; Balbus 1988). We place the local frame at the galactocentric distance R and let it corotate with the arm at the pattern speed ${{\rm{\Omega }}}_{{\rm{p}}}$. Let i denote the pitch angle of the arm segment. By taking the local approximation ($| x| ,| y| \ll R$) and assuming that the arm is tightly wound ($\mathrm{sin}i\ll 1$), the background gas velocity in the local frame arising from galaxy rotation is given by

$$\begin{eqnarray}{\boldsymbol{v}}{z}_{c}\equiv {u}_{c}\hat{{\boldsymbol{x}}}+{v}_{c}\hat{{\boldsymbol{y}}} & = & R({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{p}}})\mathrm{sin}i\;\hat{{\boldsymbol{x}}}\\ & & +[R({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{p}}})-{q}_{c}{\rm{\Omega }}x]\;\hat{{\boldsymbol{y}}},\end{eqnarray} \tag{ 1 }$$

where Ω is the rotational angular velocity at R and ${q}_{c}\equiv -d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}R$ is the local shear rate in the absence of the arm: the term involving q_c handles differential rotation in our local models. The corresponding epicycle frequency is $\kappa ={({R}^{-3}d[{{\rm{\Omega }}}^{2}{R}^{4}]/{dR})}^{1/2}={(4-2{q}_{c})}^{1/2}{\rm{\Omega }}$. Note that ${{\boldsymbol{v}}}_{c}$ in Equation (1) depends only on x and is independent of y (see also Kim & Ostriker 2002), which allows to explore the behaviors of periodic waves in the y-direction.

Since the length scales involved in gas dynamics over galactic disks are enormously large compared to electrical diffusion scales (Roberts & Yuan 1970; see also Shu 1992), the magnetic Reynolds number is much larger than unity. In this case, one can make the "frozen-in-field" approximation in which the ISM is assumed to be tightly coupled to the magnetic fields. Expanding the compressible, ideal MHD equations in the local frame, one obtains

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{{\boldsymbol{v}}}_{T})=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{{\boldsymbol{v}}}_{T}\cdot {\rm{\nabla }}{\boldsymbol{v}} & = & -{c}_{{\rm{s}}}^{2}{\rm{\nabla }}\mathrm{ln}{\rm{\Sigma }}+{q}_{c}{\rm{\Omega }}{v}_{x}\hat{{\boldsymbol{y}}}-2{\boldsymbol{\Omega }}\times {\boldsymbol{v}}\\ & & -{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{s}}}+\displaystyle \frac{1}{4\pi {\rm{\Sigma }}}({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({{\boldsymbol{v}}}_{T}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 4 }$$

together with

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Sigma }}=\int \rho {dz}$ is the gas surface density, ${\boldsymbol{v}}$ is the gas velocity induced by the arms, and ${{\boldsymbol{v}}}_{T}={\boldsymbol{v}}+{{\boldsymbol{v}}}_{c}$ is the total gas velocity in the local frame. In Equations (3) and (4), ${\boldsymbol{B}}={H}^{1/2}{{\boldsymbol{B}}}_{3{\rm{D}}}$, where ${{\boldsymbol{B}}}_{3{\rm{D}}}$ is the midplane value of the three-dimensional magnetic fields and $H={\rm{\Sigma }}/\rho$ is the disk thickness that is assumed to be constant over R (e.g., Kim & Ostriker 2001, 2002).

In Equation (3), ${{\rm{\Phi }}}_{{\rm{s}}}$ denotes the gravitational potential of the stellar spiral arms, for which we take a simple sinusoidal shape:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{s}}}={{\rm{\Phi }}}_{0}\mathrm{cos}\left(\displaystyle \frac{2\pi x}{L}\right),\end{eqnarray} \tag{ 6 }$$

constant along the y-direction. Here, ${{\rm{\Phi }}}_{0}(\lt 0)$ is the amplitude and $L=2\pi R\mathrm{sin}i/m$ is the arm-to-arm distance for m-armed spirals. We consider a domain with $-1/2\leqslant x/L\leqslant 1/2$, so that the spiral potential achieves its minimum at the domain center (x = 0). The arm strength can be characterized by the dimensionless parameter

$$\begin{eqnarray}&&{\mathcal{F}}\equiv \displaystyle \frac{m}{\mathrm{sin}i}\left(\displaystyle \frac{| {{\rm{\Phi }}}_{0}| }{{R}^{2}{{\rm{\Omega }}}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

which corresponds to the maximum gravitational force due to the spiral arms relative to the centrifugal force of the background galaxy rotation (e.g., Roberts 1969).

In the absence of spiral arms, the gaseous disk has a uniform surface density ${{\rm{\Sigma }}}_{c}$ and uniform thickness-normalized magnetic fields ${B}_{c}\hat{{\boldsymbol{y}}}$ parallel to the arms. We quantify the strength of the magnetic fields using the plasma parameter

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{v}_{{\rm{A}},c}^{2}}=\displaystyle \frac{4\pi {c}_{{\rm{s}}}^{2}{{\rm{\Sigma }}}_{c}}{{B}_{c}^{2}}=\displaystyle \frac{4\pi {P}_{c}}{{B}_{3{\rm{D}},c}^{2}},\end{eqnarray} \tag{ 8 }$$

where ${v}_{{\rm{A}},c}^{2}\equiv {B}_{c}^{2}/4\pi {{\rm{\Sigma }}}_{c}$ is the Alfvén speed, and ${P}_{c}={c}_{{\rm{s}}}^{2}{\rho }_{c}$ is the mean thermal pressure of the ISM at the galactic plane. Adopting the fiducial values ${P}_{c}/{k}_{B}\sim 2000-3000$ K cm⁻³ (Heiles 2001) and ${B}_{3{\rm{D}},c}=1.4\mu {\rm{G}}$ (Rand & Lyne 1994) in the solar neighborhood, $\beta \simeq 4$, but we consider a range of β to study various situations with differing field strength.

In Appendix A, we combine Equations (2)–(5) to obtain

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{{\boldsymbol{v}}}_{T}\cdot {\rm{\nabla }}\right){\boldsymbol{\xi }}=\displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}}{4\pi {\rm{\Sigma }}}\left(\displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}}{{\rm{\Sigma }}}\right),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{\xi }}\equiv \displaystyle \frac{{\rm{\nabla }}\times {{\boldsymbol{v}}}_{T}+2{\boldsymbol{\Omega }}}{{\rm{\Sigma }}}\end{eqnarray} \tag{ 10 }$$

is PV. This states that PV is, in general, not conserved in the presence of magnetic fields. This is unlike in 2D unmagnetized flows where ${\boldsymbol{\xi }}$ remains constant along a given streamline.

Equations (2)–(8) can be completely specified by seven parameters: q_c, m, $\mathrm{sin}i$, ${{\rm{\Omega }}}_{{\rm{p}}}/{\rm{\Omega }}$, ${c}_{{\rm{s}}}/(R{\rm{\Omega }})$, ${\mathcal{F}}$, and β. We fix to q_c = 1, m = 2, $\mathrm{sin}i=0.1$, ${{\rm{\Omega }}}_{{\rm{p}}}/{\rm{\Omega }}=0.5$, and ${c}_{{\rm{s}}}/(R{\rm{\Omega }})=0.027$, and vary ${\mathcal{F}}=5\%$–10% and $\beta =1-\infty$ for our presentation below.

2.2. Equilibrium Shock Profiles

We now seek for 1D steady solutions of Equations (2)–(5) subject to a spiral-arm forcing with strength ${\mathcal{F}}$, which will be used as an unperturbed equilibrium state in Section 3. Let us denote the steady solutions using the subscript "0" as ${{\rm{\Sigma }}}_{0}(x)$, ${{\boldsymbol{v}}}_{0}={u}_{0}(x)\hat{{\boldsymbol{x}}}+{v}_{0}(x)\hat{{\boldsymbol{y}}}$, and ${{\boldsymbol{B}}}_{0}={B}_{0}(x)\hat{{\boldsymbol{y}}}$. Then, the steady-state conditions yield

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}{u}_{T0}={{\rm{\Sigma }}}_{c}{u}_{c}=\mathrm{constant},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{u}_{T0}\displaystyle \frac{{{du}}_{0}}{{dx}}=-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{{\rm{\Sigma }}}_{0}}\displaystyle \frac{d{{\rm{\Sigma }}}_{0}}{{dx}}+2{\rm{\Omega }}{v}_{0}-\displaystyle \frac{d{{\rm{\Phi }}}_{{\rm{s}}}}{{dx}}-\displaystyle \frac{{B}_{0}}{4\pi {{\rm{\Sigma }}}_{0}}\displaystyle \frac{{{dB}}_{0}}{{dx}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{u}_{T0}\displaystyle \frac{{{dv}}_{0}}{{dx}}=-\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}{u}_{0},\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{B}_{0}{u}_{T0}={B}_{c}{u}_{c}=\mathrm{constant}.\end{eqnarray} \tag{ 14 }$$

Equations (11) and (14) imply that ${B}_{0}\propto {{\rm{\Sigma }}}_{0}$, a usual relation resulting from the conservation of mass and magnetic flux in 1D. Eliminating ${{\rm{\Sigma }}}_{0}$ and B₀ in favor of ${u}_{T0}$, Equation (12) for m = 2 reduces to

$$\begin{eqnarray}&&\left({u}_{T0}-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{u}_{T0}}-\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)\displaystyle \frac{{{du}}_{T0}}{{dx}}=2{\rm{\Omega }}{v}_{0}+R{{\rm{\Omega }}}^{2}{\mathcal{F}}\mathrm{sin}\left(\displaystyle \frac{2\pi x}{L}\right),\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{v}_{{\rm{A}}0}\equiv \displaystyle \frac{{B}_{0}}{\sqrt{4\pi {{\rm{\Sigma }}}_{0}}}={v}_{{\rm{A}},c}{\left(\displaystyle \frac{{u}_{c}}{{u}_{T0}}\right)}^{1/2},\end{eqnarray} \tag{ 16 }$$

is the Alfvén speed in an equilibrium configuration. Note that ${v}_{{\rm{A}}0}\propto {B}_{0}^{1/2}\propto {u}_{T0}^{-1/2}$.

It can be shown that the PV in the equilibrium flows is

$$\begin{eqnarray}&&{\xi }_{0}=\displaystyle \frac{| {\rm{\nabla }}\times {{\boldsymbol{v}}}_{T0}+2{\boldsymbol{\Omega }}| }{{{\rm{\Sigma }}}_{0}}=\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}{{\rm{\Sigma }}}_{c}},\end{eqnarray} \tag{ 17 }$$

which is constant everywhere even in magnetized spiral shocks. This results from the fact that an equilibrium shock satisfies ${\rm{\nabla }}\times {{\boldsymbol{B}}}_{0}=0$ (Equation (9)). Note that Equation (17) requires that the local shear rate in a steady-state spiral shock should vary as

$$\begin{eqnarray}&&{q}_{0}\equiv -\displaystyle \frac{1}{{\rm{\Omega }}}\displaystyle \frac{{{dv}}_{T0}}{{dx}}=2-(2-{q}_{c})\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{{{\rm{\Sigma }}}_{c}},\end{eqnarray} \tag{ 18 }$$

indicating that shear is reversed wherever ${{\rm{\Sigma }}}_{0}/{{\rm{\Sigma }}}_{c}\geqslant 2/(2-{q}_{c})=2$ for q_c = 1 (Balbus & Cowie 1985; Kim & Ostriker 2002).

The equilibrium velocity profiles can be obtained by solving Equations (13) and (15) simultaneously. We follow the method given by Shu et al. (1973) (see also Lee & Shu 2012; Paper I). The detailed procedure is provided in Appendix B. Figure 1 plots exemplary profiles of equilibrium spiral shocks for ${\mathcal{F}}=5\%$ and 10% and $\beta =\infty ,10,3$, and 1. A filled circle marks the magnetosonic point, ${x}_{\mathrm{mp}}$, in each shock profile. Table 1 lists the associated values of ${x}_{\mathrm{mp}}$, the shock position ${x}_{\mathrm{sh}}$, and the preshock and postshock surface densities ${{\rm{\Sigma }}}_{0}^{s-}$ and ${{\rm{\Sigma }}}_{0}^{s+}$, the density compression factor defined by

$$\begin{eqnarray}&&\mu \equiv \displaystyle \frac{{{\rm{\Sigma }}}_{0}^{s+}}{{{\rm{\Sigma }}}_{0}^{s-}}=\displaystyle \frac{{u}_{T0}^{s-}}{{u}_{T0}^{s+}},\end{eqnarray} \tag{ 19 }$$

as well as the preshock Mach number ${\mathcal{M}}\equiv {u}_{T0}^{s-}/{c}_{{\rm{s}}}$. It is apparent that magnetic fields make the shock weaker by providing magnetic pressure, reducing μ considerably. Since ${{\rm{\Sigma }}}_{0}^{s-}$ and ${\mathcal{M}}$ are insensitive to β for fixed ${\mathcal{F}}$, the reduction of the compression factor due to magnetic fields occurs primarily by making the shock front move toward the upstream direction. Note that shear is reversed in the regions behind the shock front where ${{\rm{\Sigma }}}_{0}/{{\rm{\Sigma }}}_{c}\geqslant 2$. The degree of shear reversal defined as q₀ in Equation (18) is larger as the shock becomes stronger, which makes the structures that develop as a consequence of the WI more perpendicular to the arms, as will be shown in Section 4.3.

Zoom In Zoom Out Reset image size

Table 1.  Properties of Equilibrium Spiral Shocks

${\mathcal{F}}$	β	${x}_{\mathrm{mp}}/L$	${x}_{\mathrm{sh}}/L$	${{\rm{\Sigma }}}_{0}^{s-}/{{\rm{\Sigma }}}_{c}$	${{\rm{\Sigma }}}_{0}^{s+}/{{\rm{\Sigma }}}_{c}$	μ	${\mathcal{M}}$
	$\infty$	0.007	−0.069	0.546	6.31	11.6	3.40
	100	0.008	−0.070	0.545	6.12	11.2	3.41
0.05	10	0.017	−0.078	0.538	5.01	9.29	3.45
	5	0.030	−0.085	0.534	4.34	8.12	3.48
	3	0.057	−0.093	0.528	3.79	7.17	3.52
	1	0.193	−0.125	0.512	2.62	5.12	3.63
	$\infty$	0.115	−0.005	0.316	10.9	34.5	5.87
	100	0.118	−0.006	0.316	10.4	32.7	5.88
0.10	10	0.143	−0.013	0.314	7.80	24.8	5.92
	5	0.165	−0.019	0.312	6.54	20.9	5.94
	3	0.187	−0.026	0.311	5.58	17.9	5.97
	1	0.249	−0.054	0.310	3.70	11.9	5.99

Note. For the arm and galaxy parameters of q_c = 1, m = 2, $\mathrm{sin}i=0.1$, ${{\rm{\Omega }}}_{{\rm{p}}}/{\rm{\Omega }}=0.5$, and ${c}_{{\rm{s}}}/(R{\rm{\Omega }})=0.027$. The spiral potential is minimized at x = 0.

Download table as:  ASCIITypeset image

3.1. Perturbation Equation

We apply small-amplitude Eulerian perturbations, denoted by ${{\rm{\Sigma }}}_{1}$, u₁, v₁, and ${\boldsymbol{B}}$₁, to an equilibrium shock found in the preceding section. We then linearize Equations (2)–(5) to obtain

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}}{{Dt}}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{1}}{{{\rm{\Sigma }}}_{0}}\right)+\displaystyle \frac{\partial {u}_{1}}{\partial x}+\displaystyle \frac{d\mathrm{ln}{{\rm{\Sigma }}}_{0}}{{dx}}{u}_{1}+\displaystyle \frac{\partial {v}_{1}}{\partial y}=0,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}{u}_{1}}{{Dt}}+\displaystyle \frac{{{du}}_{0}}{{dx}}{u}_{1}+{c}_{{\rm{s}}}^{2}\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{1}}{{{\rm{\Sigma }}}_{0}}\right)\\ &&-2{\rm{\Omega }}{v}_{1}-\displaystyle \frac{{B}_{0}}{4\pi {{\rm{\Sigma }}}_{0}}\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}\right){m}_{1}\\ &&-\displaystyle \frac{1}{4\pi {{\rm{\Sigma }}}_{0}}\displaystyle \frac{{{dB}}_{0}}{{dx}}\left(\displaystyle \frac{\partial {m}_{1}}{\partial x}+{B}_{0}\displaystyle \frac{{{\rm{\Sigma }}}_{1}}{{{\rm{\Sigma }}}_{0}}\right)=0,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}{v}_{1}}{{Dt}}+{c}_{{\rm{s}}}^{2}\displaystyle \frac{\partial }{\partial y}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{1}}{{{\rm{\Sigma }}}_{0}}\right)+\left(\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\right)\displaystyle \frac{{u}_{c}}{{u}_{T0}}{u}_{1}-\displaystyle \frac{1}{4\pi {{\rm{\Sigma }}}_{0}}\displaystyle \frac{{{dB}}_{0}}{{dx}}\displaystyle \frac{\partial {m}_{1}}{\partial y}=0,\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}{m}_{1}}{{Dt}}={B}_{0}{u}_{1},\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}}{{Dt}}=\displaystyle \frac{\partial }{\partial t}+{u}_{T0}\displaystyle \frac{\partial }{\partial x}+{v}_{T0}\displaystyle \frac{\partial }{\partial y},\end{eqnarray} \tag{ 24 }$$

and ${m}_{1}(x,y,t)$ is the perturbed vector potential defined through ${{\boldsymbol{B}}}_{1}\equiv {\rm{\nabla }}\times ({m}_{1}\hat{{\boldsymbol{z}}})$.

Since the coefficients in Equations (20)–(23) depend only on x, we may consider perturbations of the form

$$\begin{eqnarray}&&\left(\begin{array}{c}{{\rm{\Sigma }}}_{1}/{{\rm{\Sigma }}}_{0}\\ {u}_{1}\\ {v}_{1}\\ {m}_{1}/{B}_{c}\end{array}\right)=\left(\begin{array}{c}{S}_{1}(x)\\ {U}_{1}(x)\\ {V}_{1}(x)\\ {M}_{1}(x)\end{array}\right)\mathrm{exp}(-i\omega t+{{ik}}_{y}y),\end{eqnarray} \tag{ 25 }$$

where ω and k_y are the frequency and y-wavenumber of the perturbations, respectively. Equations (20)–(23) then reduce to

$$\begin{eqnarray}({u}_{T0}^{2}-{c}_{{\rm{s}}}^{2}-{v}_{{\rm{A}}0}^{2})\displaystyle \frac{{{dS}}_{1}}{{dx}} & = & \left[i{\omega }_{D}\left({u}_{T0}-\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)+\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\displaystyle \frac{{{du}}_{0}}{{dx}}\right]{S}_{1}\\ & & +\left(2\displaystyle \frac{{{du}}_{0}}{{dx}}-i{\omega }_{D}\right)\left(1+\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}^{2}}\right){U}_{1}\\ & & -\left[{{ik}}_{y}\left({u}_{T0}-\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)+2{\rm{\Omega }}\right]{V}_{1}\\ & & +\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{c}{u}_{T0}}\left({\omega }_{D}^{2}+2i{\omega }_{D}\displaystyle \frac{{{du}}_{0}}{{dx}}\right.\\ & & +\left.{{ik}}_{y}{u}_{T0}\displaystyle \frac{{{dv}}_{T0}}{{dx}}+{k}_{y}^{2}{u}_{T0}^{2}\right){M}_{1},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}({u}_{T0}^{2}-{c}_{{\rm{s}}}^{2}-{v}_{{\rm{A}}0}^{2})\displaystyle \frac{{{dU}}_{1}}{{dx}} & = & -\left({{ic}}_{{\rm{s}}}^{2}{\omega }_{D}+{v}_{{\rm{A}}0}^{2}\displaystyle \frac{{{du}}_{0}}{{dx}}\right){S}_{1}\\ & & +\left({{ik}}_{y}{c}_{{\rm{s}}}^{2}+2{u}_{T0}{\rm{\Omega }}\right){V}_{1}\\ & & +\left[i{\omega }_{D}\left({u}_{T0}+\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)\right.\\ & & -\left.\left({u}_{T0}+\displaystyle \frac{{c}_{{\rm{s}}}^{2}+3{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)\displaystyle \frac{{{du}}_{0}}{{dx}}\right]{U}_{1}\\ & & -\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{c}}\left({\omega }_{D}^{2}+2i{\omega }_{D}\displaystyle \frac{{{du}}_{0}}{{dx}}\right.\\ & & +\left.{{ik}}_{y}{u}_{T0}\displaystyle \frac{{{dv}}_{T0}}{{dx}}+{k}_{y}^{2}{u}_{T0}^{2}\right){M}_{1},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}{u}_{T0}\displaystyle \frac{{{dV}}_{1}}{{dx}} & = & -{{ik}}_{y}{c}_{{\rm{s}}}^{2}{S}_{1}-\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\displaystyle \frac{{u}_{c}}{{u}_{T0}}{U}_{1}\\ & & +i{\omega }_{D}{V}_{1}-{{ik}}_{y}\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{c}}\displaystyle \frac{{{du}}_{0}}{{dx}}{M}_{1},\end{eqnarray} \tag{ 28 }$$

and

$$\begin{eqnarray}&&{u}_{T0}\displaystyle \frac{{{dM}}_{1}}{{dx}}=i{\omega }_{D}{M}_{1}+\displaystyle \frac{{u}_{c}}{{u}_{T0}}{U}_{1},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{\omega }_{D}(x)=\omega -{k}_{y}{v}_{T0}\end{eqnarray} \tag{ 30 }$$

is the Doppler-shifted frequency. We take a convention that ${k}_{y}$ is a pure real number and ω is a complex number.

3.2. Shock Jump Conditions

The applied perturbations also disturb the shock front into a sinusoidal shape. Let the shape of the perturbed shock front be described by

$$\begin{eqnarray}&&{\zeta }_{1}(x,y,t)={Z}_{1}\mathrm{exp}(-i\omega t+{{ik}}_{y}y),\end{eqnarray} \tag{ 31 }$$

with the constant amplitude Z₁. The unit vectors normal and tangential to the instantaneous shock front are given by $\hat{{\boldsymbol{n}}}=(1,-{{ik}}_{y}{\zeta }_{1})$ and $\hat{{\boldsymbol{t}}}=({{ik}}_{y}{\zeta }_{1},1)$, respectively, while the velocity of the shock front is ${{\boldsymbol{v}}}_{\mathrm{sh}}=(-i\omega {\zeta }_{1},0)$ to the first order in ${\zeta }_{1}$ (Dwarkadas & Balbus 1996; Lee & Shu 2012; Paper I). It is then straightforward to show that the perturbations at the perturbed shock positions can be written as

$$\begin{eqnarray}&&{\rm{\Sigma }}({x}_{\mathrm{sh}}+{\zeta }_{1})\approx {{\rm{\Sigma }}}_{0}+{{\rm{\Sigma }}}_{1}+{\zeta }_{1}\displaystyle \frac{d{{\rm{\Sigma }}}_{0}}{{dx}},\end{eqnarray} \tag{ 32a }$$

$$\begin{eqnarray}&&{v}_{\perp }({x}_{\mathrm{sh}}+{\zeta }_{1})\approx {u}_{T0}+{u}_{1}+{\zeta }_{1}\displaystyle \frac{{{du}}_{T0}}{{dx}}+i{\omega }_{D}{\zeta }_{1},\end{eqnarray} \tag{ 32b }$$

$$\begin{eqnarray}&&{v}_{\parallel }({x}_{\mathrm{sh}}+{\zeta }_{1})\approx {v}_{T0}+{v}_{1}+{\zeta }_{1}\displaystyle \frac{{{dv}}_{T0}}{{dx}}+{{ik}}_{y}{\zeta }_{1}{u}_{T0},\end{eqnarray} \tag{ 32c }$$

$$\begin{eqnarray}&&{B}_{\perp }({x}_{\mathrm{sh}}+{\zeta }_{1})\approx {{ik}}_{y}{m}_{1}-{{ik}}_{y}{\zeta }_{1}{B}_{0},\end{eqnarray} \tag{ 32d }$$

$$\begin{eqnarray}&&{B}_{\parallel }({x}_{\mathrm{sh}}+{\zeta }_{1})\approx {B}_{0}-\displaystyle \frac{\partial {m}_{1}}{\partial x}+{\zeta }_{1}\displaystyle \frac{{{dB}}_{0}}{{dx}},\end{eqnarray} \tag{ 32e }$$

where the "⊥" and "∥" signs denote the components perpendicular and parallel to the instantaneous shock front in the stationary shock frame, respectively. All the quantities in the right-hand side (rhs) of Equation (32a) are evaluated at $x={x}_{\mathrm{sh}}$.

The jump conditions at the perturbed shock location are given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({v}_{\perp }{\rm{\Sigma }}\right)=0,\end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left(({c}_{{\rm{s}}}^{2}+{v}_{\perp }^{2}){\rm{\Sigma }}-\displaystyle \frac{{B}_{\perp }^{2}-{B}_{\parallel }^{2}}{8\pi }\right)=0,\end{eqnarray} \tag{ 33b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({\rm{\Sigma }}{v}_{\perp }{v}_{\parallel }-\displaystyle \frac{{B}_{\perp }{B}_{\parallel }}{4\pi }\right)=0,\end{eqnarray} \tag{ 33c }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({B}_{\perp }{v}_{\parallel }-{B}_{\parallel }{v}_{\perp }\right)=0,\end{eqnarray} \tag{ 33d }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({B}_{\perp }\right)=0,\end{eqnarray} \tag{ 33e }$$

(e.g., Shu 1992). Substituting Equation (32) into Equation (33), one can see that the zeroth-order terms results in Equation (64). The first-order terms are grouped to yield

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}{u}_{T0}{{\rm{\Delta }}}_{{\rm{s}}}\left({S}_{1}\right)+{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Sigma }}}_{0}{U}_{1}\right)+{{iZ}}_{1}{\omega }_{D}^{s}{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Sigma }}}_{0}\right)=0,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{{u}_{T0}^{2}+{c}_{{\rm{s}}}^{2}}{2{u}_{T0}}\right){{\rm{\Delta }}}_{{\rm{s}}}\left({S}_{1}\right)+\left(1-\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{2{u}_{T0}^{2}}\right){{\rm{\Delta }}}_{{\rm{s}}}\left({U}_{1}\right)-i{\omega }_{D}^{s}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{2{u}_{T0}^{2}\beta }{{\rm{\Delta }}}_{{\rm{s}}}\left({M}_{1}\right)\\ &&+{Z}_{1}{{\rm{\Delta }}}_{{\rm{s}}}\left[\left(\displaystyle \frac{{u}_{T0}^{2}-{c}_{{\rm{s}}}^{2}-{v}_{{\rm{A}}0}^{2}}{2{u}_{T0}^{2}}\right)\displaystyle \frac{{{du}}_{T0}}{{dx}}\right]=0,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({V}_{1}\right)-{Z}_{1}\left(\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\displaystyle \frac{{u}_{c}}{{c}_{{\rm{s}}}^{2}}-{{ik}}_{y}\right){{\rm{\Delta }}}_{{\rm{s}}}\left({u}_{T0}\right)+{{ik}}_{y}{Z}_{1}{{\rm{\Delta }}}_{{\rm{s}}}\left(\displaystyle \frac{{v}_{{\rm{A}}0}^{2}}{{u}_{T0}}\right)\\ &&-{{ik}}_{y}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{\beta }{{\rm{\Delta }}}_{{\rm{s}}}\left(\displaystyle \frac{{M}_{1}}{{u}_{T0}}\right)=0,\end{eqnarray} \tag{ 36 }$$

and

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({M}_{1}\right)-{Z}_{1}{{\rm{\Delta }}}_{{\rm{s}}}\left({B}_{0}\right)/{B}_{c}=0,\end{eqnarray} \tag{ 37 }$$

where ${\omega }_{D}^{s}={\omega }_{D}({x}_{\mathrm{sh}})$. Note that Equations (33d) and (33e) give the same result, Equation (37).

3.3. Expansion near the Magnetosonic Point

Since the left-hand sides of Equations (26) and (27) become identically zero at the magnetosonic point, the rhs should also vanish at $x={x}_{\mathrm{mp}}$ in order for regular solutions to exist. Let us expand the perturbation variables near $x={x}_{\mathrm{mp}}$ as

$$\begin{eqnarray}&&{S}_{1}={a}_{0}+{a}_{1}\eta +{\mathcal{O}}({\eta }^{2}),\end{eqnarray} \tag{ 38a }$$

$$\begin{eqnarray}&&{U}_{1}={b}_{0}+{b}_{1}\eta +{\mathcal{O}}({\eta }^{2}),\end{eqnarray} \tag{ 38b }$$

$$\begin{eqnarray}&&{V}_{1}={c}_{0}+{c}_{1}\eta +{\mathcal{O}}({\eta }^{2}),\end{eqnarray} \tag{ 38c }$$

$$\begin{eqnarray}&&{M}_{1}={d}_{0}+{d}_{1}\eta +{\mathcal{O}}({\eta }^{2}),\end{eqnarray} \tag{ 38d }$$

with coefficients ${a}_{\mathrm{0,1}}$, ${b}_{\mathrm{0,1}}$, ${c}_{\mathrm{0,1}}$, and ${d}_{\mathrm{0,1}}$. Substituting Equations (58) and (38) into Equations (26)–(29) and taking zeroth- and first-order terms in η, one can obtain a system of five linear equations for the coefficients. This indicates that the solutions near $x={x}_{\mathrm{mp}}$ can be completely specified by three constants a₀, b₀ and d₀. Since the resulting equations are not illuminating, we do not present them here.

As in Paper I, we solve Equations (26)–(29) as a boundary value problem with eigenvector $({S}_{1},{U}_{1},{V}_{1},{M}_{1},{Z}_{1})$ and eigenvalue ω. Since this is a linear problem, we are allowed to take the amplitude of one variable arbitrarily, for which we fix $\mathrm{Re}({a}_{0})=\mathrm{Im}({a}_{0})=1$ at the magnetosonic point. We choose three trial complex values for ω, b₀, and d₀. This specifies the values of S₁, U₁, V₁, and M₁ as well as their derivatives at $x={x}_{\mathrm{mp}}$ from Equation (38). We then integrate Equations (26)–(29) from $x={x}_{\mathrm{mp}}$ both in the downstream direction up to $x={x}_{\mathrm{sh}}+L$ and in the upstream direction to $x={x}_{\mathrm{sh}}$, and apply the periodic conditions to the perturbation variables. At the shock front, the fourth boundary condition (Equation (37)) determines Z₁, which is in turn used to check the other three boundary conditions. If Equations (34)–(36) are not satisfied within tolerance, we return to the first step and repeat the calculations by changing b₀, d₀, and ω, one by one, until all the perturbed jump conditions are met.

4.1. 1D Modes with ${k}_{y}=0$

We first apply the method described above to 1D perturbations with ${k}_{y}=0$. Tables 2 and 3 list the 10 lowest eigenfrequencies with differing β for ${\mathcal{F}}=5\%$ and ${\mathcal{F}}=10\%$, respectively. The modes are numbered in ascending order of $\mathrm{Re}(\omega )$. Strongly magnetized shocks with small β possess a small number of modes. In most cases, spiral shocks have a pure decaying mode (with $\mathrm{Re}(\omega )=0$ and $\mathrm{Im}(\omega )\lt 0$), a single overstable mode (with $\mathrm{Re}(\omega )\ne 0$ and $\mathrm{Im}(\omega )\gt 0$), and multiple underdamping modes (with $\mathrm{Im}(\omega )\lt 0$). The behavior of the growth rate of the overstable mode with β is not simple. In models with ${\mathcal{F}}=5\%$, magnetic fields tend to stabilize the 1D overstable mode, suppressing it completely when $\beta \leqslant 3$. In models with ${\mathcal{F}}=10\%$, on the other hand, $\mathrm{Im}(\omega )$ of the overstable mode is almost independent of β. Figure 2 compares the eigenfunction S₁ between $\beta =\infty$ and 10 cases for five odd-n modes, showing that magnetic fields do not much affect the shapes of the eigenfunctions. The number of nodes in S₁ is $2(n-3)$ for $n\geqslant 5$ regardless of β.

Zoom In Zoom Out Reset image size

Table 2.  Eigenfrequencies of One-dimensional Perturbations for ${\mathcal{F}}=0.05$

	$\beta =\infty$		$\beta =100$		$\beta =10$
Mode	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$
1	0.000	$-3.727\times {10}^{-1}$	0.000	$-3.746\times {10}^{-1}$	0.000	$-3.760\times {10}^{-1}$
2	0.692	$-1.782\times {10}^{-1}$	0.695	$-1.799\times {10}^{-1}$	0.714	$-1.859\times {10}^{-1}$
3	1.496	$+1.222\times {10}^{-3}$	1.496	$+1.123\times {10}^{-3}$	1.497	$+5.502\times {10}^{-4}$
4	2.628	$-1.380\times {10}^{-2}$	2.629	$-1.285\times {10}^{-2}$	2.640	$-1.260\times {10}^{-2}$
5	4.023	$-2.758\times {10}^{-2}$	4.037	$-2.147\times {10}^{-2}$	4.075	$-1.942\times {10}^{-2}$
6	5.554	$-2.543\times {10}^{-2}$	5.562	$-1.935\times {10}^{-2}$	5.632	$-1.610\times {10}^{-2}$
7	7.144	$-2.618\times {10}^{-2}$	7.149	$-1.614\times {10}^{-2}$	7.249	$-1.159\times {10}^{-2}$
8	8.759	$-3.077\times {10}^{-2}$	8.767	$-1.745\times {10}^{-2}$	8.898	$-8.692\times {10}^{-3}$
9	10.389	$-2.555\times {10}^{-2}$	10.406	$-1.234\times {10}^{-2}$	10.568	$-7.080\times {10}^{-3}$
10	12.033	$-2.250\times {10}^{-2}$	12.057	$-1.056\times {10}^{-2}$	12.248	$-6.936\times {10}^{-3}$
	$\beta =5$		$\beta =3$		$\beta =1$
Mode	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$
1	0.000	$-4.422\times {10}^{-1}$	0.000	$-6.030\times {10}^{-1}$	0.965	$-2.822\times {10}^{-1}$
2	0.737	$-1.920\times {10}^{-1}$	0.772	$-1.980\times {10}^{-1}$	1.501	$-1.595\times {10}^{-3}$
3	1.498	$+1.309\times {10}^{-4}$	1.499	$-2.745\times {10}^{-4}$	2.731	$-8.620\times {10}^{-3}$
4	2.652	$-1.211\times {10}^{-2}$	2.667	$-1.135\times {10}^{-2}$	4.363	$-7.434\times {10}^{-3}$
5	4.113	$-1.803\times {10}^{-2}$	4.160	$-1.595\times {10}^{-2}$	6.136	$-2.720\times {10}^{-3}$
6	5.701	$-1.381\times {10}^{-2}$	5.784	$-1.103\times {10}^{-2}$	...	...
7	7.348	$-8.594\times {10}^{-3}$	7.469	$-4.305\times {10}^{-3}$	...	...
8	9.032	$-6.711\times {10}^{-3}$	9.185	$-2.242\times {10}^{-3}$	...	...
9	10.727	$-3.669\times {10}^{-3}$	10.920	$-1.430\times {10}^{-3}$	...	...
10	12.437	$-2.614\times {10}^{-3}$	...	...	...	...

Download table as:  ASCIITypeset image

Table 3.  Eigenfrequencies of One-dimensional Perturbations for ${\mathcal{F}}=0.10$

	$\beta =\infty$		$\beta =100$		$\beta =10$
Mode	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$
1	0.354	$-6.363\times {10}^{-1}$	0.363	$-7.086\times {10}^{-1}$	0.379	$-6.983\times {10}^{-1}$
2	0.809	$-3.897\times {10}^{-1}$	0.817	$-3.963\times {10}^{-1}$	0.855	$-4.454\times {10}^{-1}$
3	1.614	$-1.078\times {10}^{-2}$	1.617	$-1.173\times {10}^{-3}$	1.619	$-1.517\times {10}^{-2}$
4	2.740	$+7.693\times {10}^{-3}$	2.741	$+8.745\times {10}^{-3}$	2.759	$+8.427\times {10}^{-3}$
5	4.316	$-3.565\times {10}^{-2}$	4.319	$-2.545\times {10}^{-2}$	4.377	$-2.297\times {10}^{-2}$
6	6.019	$-5.747\times {10}^{-2}$	6.052	$-3.498\times {10}^{-2}$	6.155	$-3.220\times {10}^{-2}$
7	7.846	$-5.831\times {10}^{-2}$	7.853	$-3.778\times {10}^{-2}$	7.996	$-2.700\times {10}^{-2}$
8	9.663	$-8.771\times {10}^{-2}$	9.678	$-3.491\times {10}^{-2}$	9.866	$-2.868\times {10}^{-2}$
9	11.509	$-7.421\times {10}^{-2}$	11.522	$-3.069\times {10}^{-2}$	11.752	$-2.190\times {10}^{-2}$
10	13.352	$-8.511\times {10}^{-2}$	13.375	$-2.640\times {10}^{-2}$	13.645	$-1.699\times {10}^{-2}$
	$\beta =5$		$\beta =3$		$\beta =1$
Mode	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$	$\mathrm{Re}(\omega )/{\rm{\Omega }}$	$\mathrm{Im}(\omega )/{\rm{\Omega }}$
1	0.894	$-5.006\times {10}^{-1}$	0.935	$-5.721\times {10}^{-1}$	1.033	$-8.362\times {10}^{-1}$
2	1.622	$-2.052\times {10}^{-2}$	1.624	$-2.678\times {10}^{-2}$	1.629	$-5.097\times {10}^{-2}$
3	2.776	$+8.532\times {10}^{-3}$	2.796	$+8.896\times {10}^{-3}$	2.877	$+9.385\times {10}^{-3}$
4	4.433	$-2.016\times {10}^{-2}$	4.498	$-1.702\times {10}^{-2}$	4.745	$-7.109\times {10}^{-3}$
5	6.250	$-2.726\times {10}^{-2}$	6.361	$-2.199\times {10}^{-2}$	6.776	$-7.802\times {10}^{-3}$
6	8.132	$-2.718\times {10}^{-2}$	8.286	$-2.087\times {10}^{-2}$	8.863	$-2.309\times {10}^{-3}$
7	10.042	$-1.719\times {10}^{-2}$	10.241	$-1.584\times {10}^{-2}$	...	...
8	11.964	$-1.949\times {10}^{-2}$	12.207	$-7.872\times {10}^{-3}$	...	...
9	13.896	$-7.560\times {10}^{-3}$	14.180	$-2.576\times {10}^{-3}$	...	...
10	...	...	...	...	...	...

Download table as:  ASCIITypeset image

To verify the growth rates of the overstable modes found by our stability analysis, we run direct MHD simulations by utilizing the Athena code (Stone et al. 2008; Stone & Gardiner 2009). Among various schemes implemented in Athena, we use the constrained corner transport method for directionally unsplit integration, the HLLE Riemann solver for flux computation (Harten et al. 1983; Einfeldt et al. 1991), and the piecewise linear method for spatial reconstruction. The simulation domain is a 1D box with length L resolved by 2048 zones. We start from an initially uniform surface density ${{\rm{\Sigma }}}_{c}$ and uniform magnetic fields with $\beta =1$ in the local frame described in Section 2. We slowly turn on the spiral arm potential amplitude and make it attain full strength ${\mathcal{F}}=10\%$ at $t=40/{\rm{\Omega }}$.

Figure 3(a) plots the temporal evolution of the gas surface density at $x/L=0$. Note that Σ increases with time as the arm strength increases and saturates at about $t=50/{\rm{\Omega }}$, after which Σ exhibits a driven-oscillator behavior. The red solid lines enveloping the density fluctuations after saturation correspond to the growth rate of $1.1\times {10}^{-2}{\rm{\Omega }}$, very close to the value given in Table 3. The inset clearly shows the oscillations of Σ over the time interval $220\leqslant t{\rm{\Omega }}\leqslant 260$. Figure 3(b) plots the Fourier-transformed power spectrum over $t{\rm{\Omega }}=200$–400. It is peaked at some specific frequencies, indicated by red arrows, equal to the real parts of the eigenfrequencies listed in Table 3. Note the dominance of the n = 3 mode with $\mathrm{Re}(\omega )/{\rm{\Omega }}=2.88$ in the power spectrum. This confirms that the results of our linear stability analysis are reliable at least for the 1D perturbations.

Zoom In Zoom Out Reset image size

It is uncertain what causes the 1D overstability of spiral shocks. It appears that the overstability arises due to a complicated interplay among various involved agents such as spiral forcing, epicycle motions, thermal pressure, magnetic fields, non-uniform background density and shear, etc. When a spiral shock is displaced from its equilibrium position, it is forced to move backward, but with an increased amplitude. Notwithstanding its nature, the growth time of the overstability amounts to $\sim 4\;\mathrm{Gyr}\;{({\rm{\Omega }}/26\;\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{kpc}}^{-1})}^{-1}$, which is much longer than that of the 2D WI presented below. This is also longer than the expected lifetime (shorter than $\sim 1\;\mathrm{Gyr}$) of spiral arms (e.g., Sellwood & Carlberg 1984; Oh et al. 2008, 2015; Speights & Westpfahl 2011, 2012). Therefore, spiral shocks can be considered stable to 1D perturbations for all practical purposes.

4.2. 2D Modes ${k}_{y}\ne 0$

4.2.1. Dispersion Relations

We now consider 2D perturbations with ${k}_{y}\ne 0$ to explore the WI of magnetized spiral shocks. Figure 4 plots the dispersion relations of ten lowest-frequency eigenmodes over $0\leqslant {k}_{y}L\leqslant 60$ for a spiral shock with ${\mathcal{F}}=5\%$ and $\beta =10$. Similarly to the 1D modes, these 2D modes are numbered in the ascending order of $\mathrm{Re}(\omega )$ at ${k}_{y}=0$. The solid and dashed lines correspond to $\mathrm{Im}(\omega )$ and $\mathrm{Re}(\omega )$, respectively. Note that $\mathrm{Re}(\omega )$ varies almost linearly with ${k}_{y},$ which indicates that the eigenmodes can be expressed as a linear superposition of entropy-vortex waves and MHD waves (Paper I). Note also that each mode becomes unstable (i.e., $\mathrm{Im}(\omega )\gt 0$) in a few ranges of ${k}_{y},$ although the corresponding growth rates for all n are less than $0.5{\rm{\Omega }}$. This is in contrast to the unmagnetized cases in which $\mathrm{Im}(\omega )$ of the n = 7 mode keeps increasing with ${k}_{y}$ (see Figure 4 of Paper I).

Zoom In Zoom Out Reset image size

To evaluate the quantitative effects of the magnetic fields on the WI, we compare the dispersion relations of overstable modes with differing β for the ${\mathcal{F}}=5\%$ and ${\mathcal{F}}=10\%$ cases in Figures 5 and 6, respectively. The growth rates and wavelengths of the most unstable modes depend on ${\mathcal{F}}$ as well as β quite sensitively. The WI of unmagnetized arms is dominated by a single (n = 7 for ${\mathcal{F}}=5\%$ and n = 4 for ${\mathcal{F}}=10\%$) mode, and this holds as long as magnetic fields are relatively weak with $\beta \geqslant 100$ for ${\mathcal{F}}=5\%$ and $\beta \geqslant 5$ for ${\mathcal{F}}=10\%$. Clearly, magnetic fields reduce both the maximum growth rate $\mathrm{Im}{(\omega )}_{\mathrm{max}}$ and the corresponding wavenumber ${k}_{y,\mathrm{max}}$, such that $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}=1.36$ and 1.13 occurring at ${k}_{y,\mathrm{max}}L=100.6$ and 78.1 for $\beta =\infty$ and 100, respectively, when ${\mathcal{F}}=5\%$. For ${\mathcal{F}}=10\%$, these values are increased to $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}=4.71$, 2.53, 1.46, and 0.90 occurring at ${k}_{y,\mathrm{max}}L=198.3$, 90.5, 47.1, and 31.9 for $\beta =\infty$, 100, 10, and 5, respectively. For more strongly magnetized arms with $\beta \leqslant 10$ for ${\mathcal{F}}=5\%$ and $\beta \leqslant 3$ for ${\mathcal{F}}=10\%$, on the other hand, there is no single dominant mode, but spiral shocks are unstable to several different modes with similar growth rates that become smaller with decreasing β and ${\mathcal{F}}$. In this case, the range of the most unstable wavenumbers is ${k}_{y,\mathrm{max}}L\sim 20$–50 for ${\mathcal{F}}=5\%$ and ${k}_{y,\mathrm{max}}L\sim 5$–30 for ${\mathcal{F}}=10\%$, largely independent of β. For arms with $\beta =1$, $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}\leqslant 0.25$ and 0.20 for ${\mathcal{F}}=5\%$ and 10%, respectively. The reduction of the growth rates is larger at larger ${k}_{y},$ suggesting that magnetic fields suppress the growth of the WI, especially for small-scale modes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 plots the eigenfunctions S₁, U₁, V₁, M₁, and ${{\rm{\Xi }}}_{1}$ of the most unstable modes with $\omega /{\rm{\Omega }}=84.4+1.13i$ and ${k}_{y}L=78.1$ for $\beta =100$ in the left panels, and with $\omega /{\rm{\Omega }}=59.1+0.50i$ and ${k}_{y}L=53.4$ for $\beta =10$ in the right panels. Also plotted as the black solid lines in the bottom panels are the solutions of Equation (54), which are in good agreement with the direct numerical results. The spiral forcing is fixed to ${\mathcal{F}}=5\%$ for both cases. When $\beta =100$, the amplitudes of the eigenfunctions decrease monotonically with the distance downstream from the shock front, similarly to the unmagnetized cases (Paper I). For more strongly magnetized spiral shocks, however, the eigenfunctions do not decay monotonically. Although $| {{\rm{\Xi }}}_{1}|$ decreases with x right after the shock front, the magnetic tension and pressure forces from the perturbed fields cause the fluctuations in the amplitude of ${{\rm{\Xi }}}_{1}$ further downstream.

Zoom In Zoom Out Reset image size

Figure 8 overlays the configuration of the perturbed magnetic fields shown as black solid lines over the real perturbed PV constructed as

$$\begin{eqnarray}&&\mathrm{Re}({\xi }_{1})=\mathrm{Re}({{\rm{\Xi }}}_{1})\mathrm{cos}({k}_{y}y)-\mathrm{Im}({{\rm{\Xi }}}_{1})\mathrm{sin}({k}_{y}y),\end{eqnarray} \tag{ 39 }$$

in the regions with $-0.2\leqslant x/L\leqslant 0.3$ and $0\leqslant y/{\lambda }_{y}\leqslant 1$, with ${\lambda }_{y}=2\pi /{k}_{y}$, for the models shown in Figure 7.³ White dots in both panels trace the wavefronts of the perturbed PV, clearly showing a discontinuity in ${k}_{x}$ across the shock located at $x/L=-0.07$ and −0.08 for $\beta =100$ and 10, respectively. Gas motions associated with entropy-vortex modes and MHD modes not only bend magnetic fields but also compress them in the postshock flows. This in turn breaks the conservation of the PV and results in a non-monotonic behavior of $| {{\rm{\Xi }}}_{1}|$. Consequently, the net reduction of the perturbed PV from one shock to next is smaller for shocks with stronger magnetic fields, which should be matched by the jump of the perturbed PV, ${{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)$, across the shock front. Note that the perturbed fields are strong inside the regions of positive $\mathrm{Re}({\xi }_{1})$ created by the WI, and reverse the directions in the regions between them.

Zoom In Zoom Out Reset image size

4.2.2. Effects of Magnetic Fields on the WI

For unmagnetized shocks, Paper I derived an analytic expression for the change of the perturbed PV across a disturbed shock, demonstrating that the shock distortion is indeed a source of the PV (see also Hayes 1957; Kevlahan 1997). Appendix C derives a similar expression by combining Equations (34)–(36) for the perturbed PV jump across a magnetized spiral shock. The derived Equation (75) recovers Equation (A8) of Paper I when $\beta =\infty$.

While ${{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)$ formally depends on all of the five perturbation variables, we find empirically that the U₁, V₁, and Z₁ terms dominate in Equation (75). Since $| {{du}}_{T0}/{dx}| \ll | {\omega }_{D}^{s}|$ for overstable modes, one can approximately write

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)}{{{ik}}_{y}} & \approx & ({{\mathcal{E}}}_{U,{\rm{H}}}+{{\mathcal{E}}}_{U,{\rm{M}}})\displaystyle \frac{{U}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}+{{\mathcal{E}}}_{V,{\rm{H}}}\displaystyle \frac{{V}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}\\ & & +({{\mathcal{E}}}_{Z,{\rm{H}}}+{{\mathcal{E}}}_{Z,{\rm{M}}})\displaystyle \frac{{Z}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}},\end{eqnarray} \tag{ 40 }$$

where the coefficients are given as

$$\begin{eqnarray}&&{{\mathcal{E}}}_{U,{\rm{H}}}=\displaystyle \frac{{(\mu -1)}^{2}}{{\mu }^{2}},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{U,{\rm{M}}}=2{\mathcal{A}}{\mathcal{B}}\left(1+\displaystyle \frac{1}{{\mu }^{2}}+\displaystyle \frac{(\mu +1){\mathcal{B}}}{\mu }\right),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{V,{\rm{H}}}=-\displaystyle \frac{{q}_{c}{\rm{\Omega }}L}{{u}_{T0}^{s-}},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{Z,{\rm{H}}}=i{\omega }_{D}^{s}\displaystyle \frac{{(\mu -1)}^{2}}{{\mu }^{2}},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{Z,{\rm{M}}}=\displaystyle \frac{2i{\omega }_{D}^{s}{\mathcal{A}}{\mathcal{B}}}{\mu (\mu +1)}\left[{(\mu +1)}^{2}-\mu (\mu -1){\mathcal{B}}\right],\end{eqnarray} \tag{ 45 }$$

with ${\mathcal{A}}$ and ${\mathcal{B}}$ defined by Equations (68) and (69), respectively. The subscripts "H" and "M" in the coefficients stand for the hydrodynamic and magnetic contributions to the perturbed PV, respectively. For unmagnetized shocks, ${\mathcal{B}}=0$, so that ${{\mathcal{E}}}_{U,{\rm{M}}}={{\mathcal{E}}}_{Z,{\rm{M}}}=0$.

Note that the U₁ term in Equation (40) results from the tangential variation of the velocity perpendicular to the shock, which plays a stabilizing role by reducing the PV by shock compression. The ${{\mathcal{E}}}_{U,{\rm{M}}}$ term explains the PV change by the perturbed magnetic pressure (Equation (33b)). That both ${{\mathcal{E}}}_{U,{\rm{H}}}$ and ${{\mathcal{E}}}_{U,{\rm{M}}}$ are always positive implies that the perturbed magnetic pressure stabilizes the WI. This can be seen more explicitly by considering a special case where the PV is retained only in U₁, with the other perturbation variables taken zero. Then, one can show that ${{\rm{\Xi }}}_{1}^{s+}/{{\rm{\Xi }}}_{1}^{s-}=(2\mu -1)/{\mu }^{2}-{{\mathcal{E}}}_{U,M}$, which is always less than unity. The V₁ term in Equation (40) comes simply from the discontinuity of the radial wavenumber across the shock, independent of the magnetic fields, which also stabilizes the WI (Paper I).

On the other hand, the Z₁ term originates from the shock deformation along the tangential direction. This is a source for the PV generation, leading to the WI. The ${{\mathcal{E}}}_{Z,{\rm{M}}}$ term is due to the perturbed magnetic stress (Equation (33c)). Since ${{\mathcal{E}}}_{Z,{\rm{H}}}$ and ${{\mathcal{E}}}_{Z,{\rm{M}}}$ always have the same sign, the stress of the deformed magnetic fields destabilizes the WI. Therefore, the role of the perturbed magnetic fields differs in the U₁ and Z₁ terms. As plotted in Figure 9, however, the ratios ${{\mathcal{E}}}_{U,{\rm{M}}}/{{\mathcal{E}}}_{U,{\rm{M}}}$ and ${{\mathcal{E}}}_{Z,{\rm{M}}}/{{\mathcal{E}}}_{Z,{\rm{H}}}$ are quite small, especially for large β and ${\mathcal{F}}$. The maximum contribution of the magnetic terms is less than 30% of the hydrodynamic terms, which occurs when ${\mathcal{F}}=5\%$ and $\beta =1$. This indicates that the effects of perturbed magnetic fields themselves to the WI are not significant.

Zoom In Zoom Out Reset image size

What is then responsible for the reduction of the growth rates in the presence of magnetic fields, as shown in Figures 5 and 6? It is through the magnetic pressure of the unperturbed background fields that tend to reduce the compression factor μ (Table 1; see also Figure 3 of Lee 2014). The amount of the PV production is smaller when a shock is weaker for fixed Z₁ (Equations (44) and (45)). More importantly, the amplitudes of the perturbation variables obtained by integrating Equations (26)–(29) decay less with x for smaller μ (e.g., Figure 7), resulting in larger values of $| {U}_{1}^{s-}|$ and $| {V}_{1}^{s-}|$ in more strongly magnetized shocks. This enhances the stabilizing role of the velocity terms relative to the destabilizing Z₁ term in Equation (40). Consequently, spiral shocks with stronger magnetic fields become more stable to the WI.

4.3. Numerical Simulation

To check the results of our 2D linear stability analysis, we run direct MHD simulations using the Athena code. As a simulation domain, we set up a rectangular box with size $L\times 2L$ that is resolved by a 2048 × 4096 grid.⁴ Initially, we construct a 1D shock profile found in Section 2.2 as a background state. We then apply small-amplitude perturbations to the background density that are realized by a flat-power Gaussian random field with a standard deviation of ${10}^{-3}{{\rm{\Sigma }}}_{0}$. We take the shearing-periodic boundary conditions at the x-boundaries and the periodic boundary conditions at the y-boundaries (e.g., Hawley et al. 1995; Kim & Ostriker 2002, 2006).

Figure 10 displays snapshots of the gas surface density in logarithmic scale and the configurations of the magnetic fields in the regions with $-0.2\leqslant x/L\leqslant 0.35$ and $0.5\leqslant y/L\leqslant 1.2$ for models with $\beta =100$, 10, 5, and 3 at $t{\rm{\Omega }}=8$, 16, 25, 40 when the WI saturates, respectively. The number of nonlinear structures along the y-direction that develop most strongly over the simulations domain is 21, 15, 12, and 9, corresponding to the wavenumber of ${k}_{y,\mathrm{max}}L=66.0$, 47.1, 37.7, and 28.3 for $\beta =100$, 10, 5, and 3, respectively. Note that the magnetic fields bend around nonlinear structures that are more strongly magnetized than the surrounding regions.

Zoom In Zoom Out Reset image size

Figure 11 compares the temporal evolution of the maximum surface density measured at $x/L=0$ in these models. The fastest growing mode in each model has a slope of 0.48, 0.22, 0.17, and 0.13, corresponding to the growth rate of $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}=1.11$, 0.50, 0.39, and 0.30, for $\beta =100$, 10, 5, and 3, respectively. These numerically measured wavelengths and growth rates are marked as star symbols in Figure 5, in good agreement with the analytic results for the n = 7, 10, 9, and 7 modes, respectively. Since various overstable modes with different n have similar maximum growth rates, which mode the system picks up should also depend on the initial power imposed by specific perturbations. In the model with $\beta =3$, for example, the initial density perturbations with ${k}_{y}L=25.1$ corresponding to the maximum growth rate of the n = 7 mode are about an order of magnitude larger than those with ${k}_{y}L=37.8$ corresponding to the peak of the n = 9 mode, emerging most strongly in the nonlinear stage, despite having a slightly ($\sim 6.7\%$) smaller growth rate.

Zoom In Zoom Out Reset image size

Finally, we remark on the level of turbulence generated by the WI. In the runs above, the density-weighted velocity dispersion, ${\sigma }_{y}={(\int {\rm{\Sigma }}{(v-{v}_{0})}^{2}{dxdy}/\int {\rm{\Sigma }}{dxdy})}^{1/2}$, in the direction parallel to the arms is found to be $\sim 1.4,1.1,0.8,$ and $0.7\;\mathrm{km}\;{{\rm{s}}}^{-1}$ at the time when the WI saturates for $\beta =100$, 10, 5, and 3, respectively, which are interestingly equal approximately to $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{k}_{y,\mathrm{max}}$.⁵ Due to nonlinear interactions and mergers of clumps created by the WI, ${\sigma }_{y}$ increases further by about a factor of 2. This suggests that the WI can be a considerable, although not dominant, source of turbulence energy in the ISM since the energy injection occurs to the densest part of the gas.

We have investigated the WI of magnetized spiral shocks in a galactic disk using both a linear stability analysis and direct MHD simulations. This is a straightforward extension of Paper I, which studied the case of unmagnetized spiral shocks. We assume that the gas disk is infinitesimally thin, isothermal, and non-self-gravitating. We parameterize the strengths of the stellar spiral arms and magnetic fields using the dimensionless parameters ${\mathcal{F}}$ and β, respectively (Equations (7) and (8)). As background states, we first obtain the steady spiral shock profiles with differing ${\mathcal{F}}$ and β. We then impose small-amplitude wiggling perturbations to the equilibrium spiral shocks, and calculate their dispersion relations as a boundary-value problem with eigenfrequencies. Our local MHD simulations readily pick up the most unstable modes for a given set up parameters, with the numerical growth rates very close to the results of the linear stability analysis.

The existence of the overstable modes proves that the WI is physical in origin, resulting from the accumulation of the perturbed PV from a distorted shock front. Our results show that magnetic fields suppress the WI, but not completely, at least for $\beta \geqslant 1$. The stabilizing role of magnetic fields is not from the perturbed fields but directly from the background unperturbed fields that tend to reduce the shock compression factor μ by exerting magnetic pressure. When magnetic fields are relatively weak, the WI is dominated by a single dominant mode. When ${\mathcal{F}}=5\%$, the most unstable mode has a growth rate $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}=1.36$ and 1.13 occurring at ${k}_{y,\mathrm{max}}L=100.6$ and 78.1 for $\beta =\infty$ and 100, respectively. When ${\mathcal{F}}=10\%$, these values are increased to $\mathrm{Im}{(\omega )}_{\mathrm{max}}/{\rm{\Omega }}=4.71$, 2.53, 1.46, and 0.90 occurring at ${k}_{y,\mathrm{max}}L=198.3$, 90.5, 47.1, and 31.9 for $\beta =\infty$, 100, 10, and 5, respectively. For more strongly magnetized arms with $\beta \leqslant 10$ for ${\mathcal{F}}=5\%$ and $\beta \leqslant 3$ for ${\mathcal{F}}=10\%$, on the other hand, several overstable modes have similar growth rates that become smaller with decreasing β, while the wavelength range of the most unstable modes is ${k}_{y,\mathrm{max}}L\sim 20$–50 for ${\mathcal{F}}=5\%$ and ${k}_{y,\mathrm{max}}L$ ∼ 5–30 for ${\mathcal{F}}=10\%$, insensitive to $\beta$.

We have found that magnetic fields stabilize the WI at small scales by reducing the shock compression factor. When $\beta =10$, the most unstable wavelength is decreased by a factor of about 2 and 4 for ${\mathcal{F}}=5\%$ and 10%, respectively, compared to the unmagnetized cases. The stabilization of the WI by magnetic fields has already been observed in previous nonlinear simulations of spiral galaxies. For instance, Kim & Ostriker (2006) ran 2D local models for the formation of feathers and provided the magnetic field topology at the nonlinear stage resulting from the vorticity generation near the shock front. They further showed that vortical clumps produced merge into massive clouds with mass $\sim {10}^{7}{M}_{\odot }$ each in the presence of self-gravity. Grid-based global simulations by Shetty & Ostriker (2006) found no indication of the development of the WI in a model with $\beta =1$ and ${\mathcal{F}}=10\%$. Based on our results, this is presumably not because magnetic fields completely suppress the WI but because its growth time is longer than the simulation time span (two orbits) in their models with very weak ($\sim 0.1\%$) initial perturbations.⁶ By running particle-based simulations, on the other hand, Dobbs & Price (2008) found non-axisymmetric structures, albeit weak, sill grow in models with $\beta =1$. This is probably because they might have large-amplitude density perturbations arising from the Poisson noises in the initial particle distributions, helping the WI readily manifest in their simulations.

It is interesting to compare the predicted wavelength ${\lambda }_{y,\mathrm{max}}=2\pi /{k}_{y,\mathrm{max}}$ from our linear stability analysis with the observed spacings of feathers. By analyzing the Hubble archival data, La Vigne et al. (2006) measured the feather spacings in grand-design spiral galaxies. They reported that M51 and M74 have the feather spacing that increases from $\sim 0.2\;\mathrm{kpc}$ in the inner regions ($R\sim 1\;\mathrm{kpc}$) to $\sim 1\;\mathrm{kpc}$ in the outer regions ($R\sim 10\;\mathrm{kpc}$). Adopting the arm pitch angles of $i=21\buildrel{\circ}\over{.} 1$ for M51 (Shetty et al. 2007) and $i=15\buildrel{\circ}\over{.} 7$ for M74 (Gusev & Efremov 2013) (see also Honig & Reid 2015), these correspond to ${\lambda }_{y,\mathrm{max}}/L\sim 0.1$–0.2, consistent with our results for magnetized arms with $\beta \lesssim 10$ for ${\mathcal{F}}=5\%$ or with $\beta \sim 5$–10 for ${\mathcal{F}}=10\%$.⁷ Elmegreen & Elmegreen (1983) presented the separations of H ii regions and H i superclouds along arms in various spiral galaxies. Taking m = 2 and $i=20^\circ$ arbitrarily, the observed separations listed in their Table 2 are distributed in the range ${\lambda }_{y,\mathrm{max}}/L=0.1$–1, with a mean value of ∼0.4, a factor of ∼ 2–4 larger than the mean feather spacings mentioned above. This may indicate that H ii regions and H i superclouds represent highly nonlinear structures created by mergers of WI-induced feathers.

Recently, Puerari et al. (2014) introduced a new method to locally determine the distributions of the pitch angles of spiral arms and their substructures. Applying the method to 8 μm Spitzer images of M51 and M81, they found that the mean difference between the pitch angles of the main spiral arms and the interarm secondary structures is ${\rm{\Delta }}i\sim 10^\circ$–30°. We apply the same method to the results of our MHD simulations to calculate the pitch angles of the nonlinear structures resulting from the WI off the main arms. Figure 12 plots ${\rm{\Delta }}i$ when the structures saturate in models with $\beta =3$–100 and ${\mathcal{F}}=5\%$, averaged over the immediate postshock regions from $x={x}_{\mathrm{sh}}$ to $x={x}_{\mathrm{sh}}+L/4$. Note that ${\rm{\Delta }}i$ is smaller in more strongly magnetized models. This is mostly because of shear reversal in the postshock regions, discussed in Section 2.2, which is larger for stronger shocks. Since shear reversal tends to decrease ${k}_{x},$ ${\rm{\Delta }}i={\mathrm{tan}}^{-1}({k}_{y}/{k}_{x})$ is smaller for weaker shocks with smaller β. The numerical values of ${\rm{\Delta }}i\sim 20^\circ$–35° for $\beta \sim 3$–10 are close to the values obtained by Puerari et al. (2014) for M51 and M81.

Zoom In Zoom Out Reset image size

As discussed above, the theoretical predictions of our analysis based on 2D, non-self-gravitating models are similar to the observed properties of feathers, suggesting that the WI may be responsible for the feather formation in the presence of magnetic fields. However, it still remains to be seen whether the WI would grow into observed interarm features in real disk galaxies for the following reasons. First, Kim & Ostriker (2006) showed that spiral shocks in vertically stratified disk exhibit non-steady flapping motions in the direction perpendicular to the arms. Together with strong vertical shear, these shock flapping motions tend to disrupt coherent vortical structures and thus prevent the growth of the WI.

Second, Lee & Shu (2012) and Lee (2014) recently conducted a parameter study of feathering instability of a magnetized, self-gravitating spiral shock. In particular, Lee (2014) found that the most unstable mode of the feathering instability is $\sim 530\;\mathrm{pc}$ in their M51 arm, similar to the observed feather spacing reported by La Vigne et al. (2006). He also found that the growth rate of the feathering instability is typically $\sim {\rm{\Omega }}$, which is larger than the growth rate (∼0.2–$0.3{\rm{\Omega }}$) of the WI for $\beta =1$ studied in this work. Without considering the effects of the background shear, however, they were unable to study the combined effects of the wiggle and feathering instabilities. It is possible that the WI grows fast and provides seeds for the onset of the feathering instability. Or, the WI is completely suppressed by shock flapping motions in stratified disks. To fully understand the feather formation, therefore, it is necessary to study the stability of spiral shocks in a vertically extended, self-gravitating, magnetized disk that harbor not only the wiggle and feathering instabilities but also the Parker instability.

We acknowledge helpful discussions with J. Kim and D. Ryu. We are also grateful to the referee for a thoughtful report. This work was supported by the National Research Foundation of Korea (NRF) grant, No. 2008-0060544, funded by the Korea government (MSIP). The computation of this work was supported by the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2014-C3-003).

Here, we derive an equation for the evolution of the perturbed PV. Equation (3) can be decomposed as

$$\begin{eqnarray}&&\displaystyle \frac{{{Du}}_{T}}{{Dt}}=-{c}_{{\rm{s}}}^{2}\displaystyle \frac{\partial \mathrm{ln}{\rm{\Sigma }}}{\partial x}+2{\rm{\Omega }}v-\displaystyle \frac{d{{\rm{\Phi }}}_{{\rm{s}}}}{{dx}}-\displaystyle \frac{{B}_{y}}{4\pi {\rm{\Sigma }}}\left(\displaystyle \frac{\partial {B}_{y}}{\partial x}-\displaystyle \frac{\partial {B}_{x}}{\partial y}\right),\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{Dv}}_{T}}{{Dt}}=-{c}_{{\rm{s}}}^{2}\displaystyle \frac{\partial \mathrm{ln}{\rm{\Sigma }}}{\partial y}-\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}-{q}_{c}{\rm{\Omega }}{u}_{T}+\displaystyle \frac{{B}_{x}}{4\pi {\rm{\Sigma }}}\left(\displaystyle \frac{\partial {B}_{y}}{\partial x}-\displaystyle \frac{\partial {B}_{x}}{\partial y}\right),\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{D}{{Dt}}=\displaystyle \frac{\partial }{\partial t}+{{\boldsymbol{v}}}_{T}\cdot {\rm{\nabla }}\end{eqnarray} \tag{ 48 }$$

is the Lagrangian time-derivative, and we have used the identity ${{\boldsymbol{v}}}_{T}\cdot {\rm{\nabla }}{\boldsymbol{v}}={{\boldsymbol{v}}}_{T}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{T}+{q}_{c}{\rm{\Omega }}{u}_{T}$.

Differentiating Equations (46) and (47) with respect to y and x, respectively, and then subtracting the resulting equations, we obtain

$$\begin{eqnarray}\displaystyle \frac{D}{{Dt}}\left({\rm{\nabla }}\times {{\boldsymbol{v}}}_{T}+2{\rm{\Omega }}\right) & = & \left({\rm{\nabla }}\times {{\boldsymbol{v}}}_{T}+2{\rm{\Omega }}\right)\displaystyle \frac{D}{{Dt}}\mathrm{ln}{\rm{\Sigma }}\\ & & +\displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}}{4\pi }\left(\displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}}{{\rm{\Sigma }}}\right),\end{eqnarray} \tag{ 49 }$$

which can be simplified to

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{\xi }}}{{Dt}}=\displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}}{4\pi {\rm{\Sigma }}}\left(\displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}}{{\rm{\Sigma }}}\right),\end{eqnarray} \tag{ 50 }$$

where ${\boldsymbol{\xi }}$ is the PV defined in Equation (10). Equation (50) states that ${\boldsymbol{\xi }}$ is in general not conserved in magnetized flows, which is in contrast to the hydrodynamic case where $D{\boldsymbol{\xi }}/{Dt}=0$.

After linearizing Equation (10), one can write the perturbed PV as

$$\begin{eqnarray}&&{\xi }_{1}=\displaystyle \frac{| {\rm{\nabla }}\times {{\boldsymbol{v}}}_{1}| }{{{\rm{\Sigma }}}_{0}}-{\xi }_{0}\displaystyle \frac{{{\rm{\Sigma }}}_{1}}{{{\rm{\Sigma }}}_{0}},\end{eqnarray} \tag{ 51 }$$

where ${\xi }_{0}$ is the PV in the background unperturbed shock (Equation (17)). Analogous to Equation (25), we define the amplitude ${{\rm{\Xi }}}_{1}$ of ${\xi }_{1}$ as

$$\begin{eqnarray}&&{{\rm{\Xi }}}_{1}\equiv {\xi }_{1}(x,y,t){e}^{i\omega t-{{ik}}_{y}y}=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{0}}\left(\displaystyle \frac{{{dV}}_{1}}{{dx}}-{{ik}}_{y}{U}_{1}\right)-{\xi }_{0}{S}_{1}.\end{eqnarray} \tag{ 52 }$$

In terms of the perturbation variables, Equation (50) then becomes

$$\begin{eqnarray}\left(-i{\omega }_{D}+{u}_{T0}\displaystyle \frac{d}{{dx}}\right){{\rm{\Xi }}}_{1} & = & \displaystyle \frac{{{ik}}_{y}{c}_{{\rm{s}}}^{2}}{\beta {{\rm{\Sigma }}}_{0}}\left[\left(\displaystyle \frac{{d}^{2}}{{{dx}}^{2}}\mathrm{ln}{{\rm{\Sigma }}}_{0}+{k}_{y}^{2}\right){M}_{1}\right.\\ & & -\left.\left(\displaystyle \frac{1}{{{\rm{\Sigma }}}_{c}}\displaystyle \frac{d{{\rm{\Sigma }}}_{0}}{{dx}}{S}_{1}+\displaystyle \frac{{d}^{2}{M}_{1}}{{{dx}}^{2}}\right)\right].\end{eqnarray} \tag{ 53 }$$

The terms in the first parentheses in the rhs of Equation (53) are due to magnetic tension, while those from magnetic pressure are given in the second parentheses.

A formal solution of Equation (53) is given by

$$\begin{eqnarray}&&{{\rm{\Xi }}}_{1}(x)=({{\rm{\Xi }}}_{1}^{s+}+{{\rm{\Psi }}}_{T}+{{\rm{\Psi }}}_{{\rm{P}}})\cdot {e}^{-\tau \mathrm{Im}(\omega )}\mathrm{exp}\left(i{\displaystyle \int }_{{x}_{\mathrm{sh}}}^{x}{k}_{x,v}(x){dx}\right)\end{eqnarray} \tag{ 54 }$$

where ${{\rm{\Xi }}}_{1}^{s+}$ is the PV at the immediate postshock regions, $\tau ={\displaystyle \int }_{{x}_{\mathrm{sh}}}^{x}{dx}/{u}_{T0}$ is the Lagrangian time elapsed from the shock front, and ${k}_{x}(x)=(\mathrm{Re}(\omega )-{v}_{T0}{k}_{y})/{u}_{T0}$ is the local x-wavenumber. In Equation (54), ${{\rm{\Psi }}}_{T}$ and ${{\rm{\Psi }}}_{{\rm{P}}}$ are defined by

$$\begin{eqnarray}{{\rm{\Psi }}}_{T} & \equiv & \displaystyle \frac{{{ik}}_{y}{c}_{{\rm{s}}}^{2}}{\beta {u}_{T0}{{\rm{\Sigma }}}_{0}}{\displaystyle \int }_{{x}_{\mathrm{sh}}}^{x}\left(\displaystyle \frac{{d}^{2}}{{dx}{^{\prime} }^{2}}\mathrm{ln}{{\rm{\Sigma }}}_{0}+{k}_{y}^{2}\right){M}_{1}\cdot {e}^{\tau (x^{\prime} )\mathrm{Im}(\omega )}\\ & & \times \mathrm{exp}\left(-i{\displaystyle \int }_{{x}_{\mathrm{sh}}}^{{x}^{\prime }}{k}_{x}(x^{\prime\prime} ){dx}^{\prime\prime} \right){dx}^{\prime} ,\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}{{\rm{\Psi }}}_{{\rm{P}}} & \equiv & -\displaystyle \frac{{{ik}}_{y}{c}_{{\rm{s}}}^{2}}{\beta {u}_{T0}{{\rm{\Sigma }}}_{0}}{\displaystyle \int }_{{x}_{\mathrm{sh}}}^{x}\left(\displaystyle \frac{1}{{{\rm{\Sigma }}}_{c}}\displaystyle \frac{d{{\rm{\Sigma }}}_{0}}{{dx}^{\prime} }{S}_{1}+\displaystyle \frac{{d}^{2}{M}_{1}}{{dx}{^{\prime} }^{2}}\right)\cdot {e}^{\tau (x^{\prime} )\mathrm{Im}(\omega )}\\ & & \times \mathrm{exp}\left(-i{\displaystyle \int }_{{x}_{\mathrm{sh}}}^{{x}^{\prime }}{k}_{x}(x^{\prime\prime} ){dx}^{\prime\prime} \right){dx}^{\prime} ,\end{eqnarray} \tag{ 56 }$$

representing the contributions to the perturbed PV due to magnetic tension and magnetic pressure, respectively.

The coefficient of ${{du}}_{T0}/{dx}$ in the left-hand side of Equation (15) vanishes at the magnetosonic point, $x={x}_{\mathrm{mp}}$, where

$$\begin{eqnarray}&&{\left({u}_{T0}-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{u}_{T0}}-\displaystyle \frac{{c}_{{\rm{s}}}^{2}{u}_{c}}{{u}_{T0}^{2}\beta }\right)}_{x={x}_{\mathrm{mp}}}=0,\end{eqnarray} \tag{ 57 }$$

which is a cubic equation in ${u}_{T0}$. A subsonic gas should accelerate to a supersonic speed through the magnetosonic point. To obtain smooth solutions, we expand ${u}_{T0}$ and v₀ around $x={x}_{\mathrm{mp}}$ as

$$\begin{eqnarray}&&{u}_{T0}={\alpha }_{0}+{\alpha }_{1}\eta +{\alpha }_{2}{\eta }^{2}+{\mathcal{O}}({\eta }^{3}),\end{eqnarray} \tag{ 58a }$$

$$\begin{eqnarray}&&{v}_{0}={\gamma }_{0}+{\gamma }_{1}\eta +{\gamma }_{2}{\eta }^{2}+{\mathcal{O}}({\eta }^{3}),\end{eqnarray} \tag{ 58b }$$

where $\eta \equiv x-{x}_{\mathrm{mp}}$ and ${\alpha }_{0}$ denotes the positive real root of Equation (57) for given u_c, ${c}_{{\rm{s}}},$ and $\beta$. The coefficients ${\alpha }_{\mathrm{1,2}}$ and ${\gamma }_{\mathrm{0,1,2}}$ are to be determined by series expansions near $\eta =0$.

Substituting Equation (58) in Equations (13) and (15) and keeping terms up to the second order in $\eta ,$ one obtains

$$\begin{eqnarray}&&{\gamma }_{0}=-\displaystyle \frac{{\mathcal{F}}}{2}R{\rm{\Omega }}\mathrm{sin}\left(\displaystyle \frac{2\pi {x}_{\mathrm{mp}}}{L}\right),\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&{\alpha }_{1}={\left(3-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\alpha }_{0}^{2}}\right)}^{-1/2}{\left[2{\rm{\Omega }}{\gamma }_{1}+\displaystyle \frac{2{\mathcal{F}}{{\rm{\Omega }}}^{2}}{\mathrm{sin}i}\mathrm{cos}\left(\displaystyle \frac{2\pi {x}_{\mathrm{mp}}}{L}\right)\right]}^{1/2},\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{\gamma }_{1}=\displaystyle \frac{{u}_{c}-{\alpha }_{0}}{{\alpha }_{0}}{\rm{\Omega }},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{\alpha }_{2}=\displaystyle \frac{1}{3}{\left(3-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\alpha }_{0}^{2}}\right)}^{-1}\left[\displaystyle \frac{{\alpha }_{1}^{2}}{{\alpha }_{0}}\left(3-2\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\alpha }_{0}^{2}}\right)-\displaystyle \frac{{u}_{c}{{\rm{\Omega }}}^{2}}{{\alpha }_{0}^{2}}+\displaystyle \frac{4{\gamma }_{0}}{{\alpha }_{1}}\displaystyle \frac{{\pi }^{2}}{{L}^{2}}{\rm{\Omega }}\right],\end{eqnarray} \tag{ 62 }$$

and

$$\begin{eqnarray}&&{\gamma }_{2}=-\displaystyle \frac{1+{\gamma }_{1}}{2{\alpha }_{0}}{\alpha }_{1}.\end{eqnarray} \tag{ 63 }$$

With these coefficients, ${u}_{T0}$ and v₀ in Equation (58) vary smoothly near the magnetosonic point.

The equilibrium spiral-shock profiles should satisfy the following jump conditions at the shock front, $x={x}_{\mathrm{sh}}$:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({u}_{T0}{{\rm{\Sigma }}}_{0}\right)=0,\end{eqnarray} \tag{ 64a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left(({c}_{{\rm{s}}}^{2}+{u}_{T0}^{2}){{\rm{\Sigma }}}_{0}+\displaystyle \frac{{B}_{0}^{2}}{8\pi }\right)=0,\end{eqnarray} \tag{ 64b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({v}_{0}\right)=0,\end{eqnarray} \tag{ 64c }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{s}}}\left({u}_{T0}{B}_{0}\right)=0,\end{eqnarray} \tag{ 64d }$$

where ${{\rm{\Delta }}}_{{\rm{s}}}(f)\equiv {f}^{s+}-{f}^{s-}$, with the superscripts "s+" and "s-" representing the quantities evaluated at the immediate behind ($x={x}_{\mathrm{sh}}+0$) and ahead ($x={x}_{\mathrm{sh}}-0$) of the shock front, respectively. Note that Equations (64a) and (64d) are satisfied automatically from Equations (11) and (14).

Equilibrium profiles of magnetized spiral shocks can be constructed as follows. For given ${\mathcal{F}}$ and β, we first choose the magnetosonic point ${x}_{\mathrm{mp}}$ arbitrarily, and integrate Equations (13) and (15) starting from $x={x}_{\mathrm{mp}}$ in both the upstream and downstream directions. We apply the periodic boundary conditions at $x/L=\pm 1/2$ and determine the shock position ${x}_{\mathrm{sh}}$ by imposing Equation (64c). We then check whether Equation (64b) is also fulfilled. If not within tolerance, we return to the first step to change ${x}_{\mathrm{mp}}$ until all the jump conditions are satisfied.

Here we derive an algebraic expression for the jump condition of the perturbed PV, ${{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)={{\rm{\Xi }}}_{1}^{s+}-{{\rm{\Xi }}}_{1}^{s-}$, at the shock front ($x={x}_{\mathrm{sh}}$) in terms of the preshock variables (${S}_{1}^{s-}$, ${U}_{1}^{s-}$, ${V}_{1}^{s-}$, and ${M}_{1}^{s-}$) and Z₁.

With the help of Equation (19), the jump conditions (Equation (64)) of the background steady shock can be combined to yield

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{u}_{T0}^{s+}}{{c}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{1}{\mu }+\displaystyle \frac{\mu +1}{2\mu \beta }\displaystyle \frac{{u}_{c}}{{u}_{T0}^{s-}},\end{eqnarray} \tag{ 65 }$$

for the postshock Mach number, and

$$\begin{eqnarray}&&{{\mathcal{M}}}^{2}={\left(\displaystyle \frac{{u}_{T0}^{s-}}{{c}_{{\rm{s}}}}\right)}^{2}=\mu +\displaystyle \frac{\mu (\mu +1)}{2\beta }\displaystyle \frac{{u}_{c}}{{u}_{T0}^{s-}},\end{eqnarray} \tag{ 66 }$$

for the preshock Mach number. For later purposes, we calculate

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{u}_{T0}^{s+}}{{dx}}-\displaystyle \frac{d\mathrm{ln}{u}_{T0}^{s-}}{{dx}}=-{\mathcal{A}}\left[{(\mu +1)}^{2}+2\mu {\mathcal{B}}\right]\displaystyle \frac{d\mathrm{ln}{u}_{T0}^{s-}}{{dx}},\end{eqnarray} \tag{ 67 }$$

where

$$\begin{eqnarray}&&{\mathcal{A}}\equiv {\left(2\mu +1-\displaystyle \frac{{\mu }^{2}}{{{\mathcal{M}}}^{2}}\right)}^{-1},\end{eqnarray} \tag{ 68 }$$

and

$$\begin{eqnarray}&&{\mathcal{B}}\equiv 1-\displaystyle \frac{\mu }{{{\mathcal{M}}}^{2}}.\end{eqnarray} \tag{ 69 }$$

Note that ${\mathcal{M}}={\mu }^{1/2}$, ${\mathcal{A}}=1/(1+\mu )$, and ${\mathcal{B}}=0$ for isothermal hydrodynamic shocks. Table 1 shows that ${\mathcal{M}}$ is insensitive to β for given ${\mathcal{F}}$. Both ${\mathcal{A}}$ and ${\mathcal{B}}$ are positive definite.

It is straightforward to show that Equation (34) becomes

$$\begin{eqnarray}&&{S}_{1}^{s+}={S}_{1}^{s-}+[{U}_{1}^{s-}-\mu {U}_{1}^{s+}-(\mu -1)i{\omega }_{D}^{s}{Z}_{1}]/{u}_{T0}^{s-},\end{eqnarray} \tag{ 70 }$$

while Equation (37) is simplified to

$$\begin{eqnarray}&&{M}_{1}^{s+}={M}_{1}^{s-}+(\mu -1)({u}_{c}/{u}_{T0}^{s-}){Z}_{1}.\end{eqnarray} \tag{ 71 }$$

On the other hand, Equation (35) utilizing Equations (65) and (71) becomes

$$\begin{eqnarray}\displaystyle \frac{{U}_{1}^{s+}}{{\mathcal{A}}} & = & -\displaystyle \frac{1}{\mu }[\mu +1+(\mu +2){\mathcal{B}}]{U}_{1}^{s-}\\ & & -\displaystyle \frac{(\mu +1){\mathcal{B}}}{\mu }\left({u}_{T0}^{s-}{S}_{1}^{s-}+\displaystyle \frac{2i{\omega }_{D}{u}_{T0}^{s-}}{{u}_{c}}{M}_{1}^{s-}\right)\\ & & -i{\omega }_{D}^{2}\left[\displaystyle \frac{{(\mu +1)}^{2}}{\mu }+(\mu -1){\mathcal{B}}\right]{Z}_{1}\\ & & +(1+\mu -{\mathcal{B}})\displaystyle \frac{\mu -1}{\mu }\displaystyle \frac{{{du}}_{T0}^{s-}}{{dx}}{Z}_{1}.\end{eqnarray} \tag{ 72 }$$

Plugging Equation (72) into Equation (70) allows to write ${S}_{1}^{s+}$ in terms of the preshock quantities. Lastly, Equation (36) together with Equation (71) gives

$$\begin{eqnarray}{V}_{1}^{s+} & = & {V}_{1}^{s-}-(\mu -1)\left[\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\displaystyle \frac{{u}_{c}}{{u}_{T0}^{s-}}-{{ik}}_{y}\left(\displaystyle \frac{{u}_{T0}^{s-}}{\mu }+\displaystyle \frac{{u}_{c}}{\beta {{\mathcal{M}}}^{2}}\right)\right]\\ & & \times {Z}_{1}+\displaystyle \frac{{{ik}}_{y}(\mu -1){c}_{{\rm{s}}}}{\beta {\mathcal{M}}}{M}_{1}^{s-}.\end{eqnarray} \tag{ 73 }$$

Combining Equations (28) and (52), the jump of the perturbed PV can be written as

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right) & = & -\displaystyle \frac{{{ik}}_{y}}{{{\rm{\Sigma }}}_{0}^{s-}}\left(\displaystyle \frac{{U}_{1}^{s+}}{\mu }-{U}_{1}^{s-}\right)\\ & & -\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\displaystyle \frac{{u}_{c}}{{({u}_{T0}^{s-})}^{2}{{\rm{\Sigma }}}_{0}^{s-}}(\mu {U}_{1}^{s+}-{U}_{1}^{s-})\\ & & -\left(\displaystyle \frac{{{ik}}_{y}{u}_{T0}^{s-}}{{{\mathcal{M}}}^{2}}+\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\displaystyle \frac{{u}_{c}}{{u}_{T0}^{s-}}\right)\left(\displaystyle \frac{{S}_{1}^{s+}}{{{\rm{\Sigma }}}_{0}^{s-}}-\displaystyle \frac{{S}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}\right)\\ & & +\displaystyle \frac{i{\omega }_{D}^{s}}{{u}_{T0}^{s-}{{\rm{\Sigma }}}_{0}^{s-}}({V}_{1}^{s+}-{V}_{1}^{s-})\\ & & -\displaystyle \frac{{{iq}}_{c}{\rm{\Omega }}{{Lk}}_{y}}{{u}_{T0}^{s-}{{\rm{\Sigma }}}_{0}^{s-}}{V}_{1}^{s-}-\displaystyle \frac{{{ik}}_{y}}{\beta {{\mathcal{M}}}^{2}{{\rm{\Sigma }}}_{0}^{s-}}\\ & & \times \left(\mu \displaystyle \frac{{{du}}_{T0}^{s+}}{{dx}}{M}_{1}^{s+}-\displaystyle \frac{{{du}}_{T0}^{s-}}{{dx}}{M}_{1}^{s-}\right).\end{eqnarray} \tag{ 74 }$$

Using Equations (70)–(73), we eliminate ${S}_{1}^{s+}$, ${U}_{1}^{s+}$, ${V}_{1}^{s+}$, and ${M}_{1}^{s+}$ in Equation (74) and arrange the terms to obtain

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)}{{{ik}}_{y}} & = & {{\mathcal{E}}}_{S}\displaystyle \frac{{S}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}+{{\mathcal{E}}}_{U}\displaystyle \frac{{U}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}+{{\mathcal{E}}}_{V}\displaystyle \frac{{V}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}\\ & & +{{\mathcal{E}}}_{{\rm{M}}}\displaystyle \frac{{M}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}}+{{\mathcal{E}}}_{Z}\displaystyle \frac{{Z}_{1}^{s-}}{{{\rm{\Sigma }}}_{0}^{s-}},\end{eqnarray} \tag{ 75 }$$

with the coefficients defined as

$$\begin{eqnarray}&&{{\mathcal{E}}}_{S}=\displaystyle \frac{(\mu +1){\mathcal{A}}{\mathcal{B}}}{{\mu }^{2}}\left(1-\displaystyle \frac{{\mu }^{2}}{{{\mathcal{M}}}^{2}}\right){u}_{T0}^{s-},\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{U}=\displaystyle \frac{{(\mu -1)}^{2}}{{\mu }^{2}}+2{\mathcal{A}}{\mathcal{B}}\left(1+\displaystyle \frac{1}{{\mu }^{2}}+\displaystyle \frac{(\mu +1){\mathcal{B}}}{\mu }\right),\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{{\mathcal{E}}}_{V}=-\displaystyle \frac{{q}_{c}{\rm{\Omega }}L}{{u}_{T0}^{s-}},\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}{{\mathcal{E}}}_{{\rm{M}}} & = & \displaystyle \frac{2{\mathcal{A}}{\mathcal{B}}}{{\mu }^{2}(1+\mu )}\displaystyle \frac{{u}_{T0}^{s-}}{{u}_{c}}\left\{\left(4{\mu }^{2}-1+\displaystyle \frac{2}{{\mathcal{A}}}\right)\right.\\ & & \times \left.\displaystyle \frac{{{du}}_{T0}^{s-}}{{dx}}+i{\omega }_{D}^{s}\left[1-{\mu }^{2}+\mu (1+\mu +2{\mu }^{2}){\mathcal{B}}\right]\right\},\end{eqnarray} \tag{ 79 }$$

and

$$\begin{eqnarray}{{\mathcal{E}}}_{Z} & = & \displaystyle \frac{{(\mu -1)}^{2}}{\mu (1+\mu )}\left(1+\displaystyle \frac{1}{{\mathcal{M}}}\right)\displaystyle \frac{{{du}}_{T0}^{s-}}{{dx}}+i{\omega }_{D}^{s}\displaystyle \frac{{(\mu -1)}^{2}}{{\mu }^{2}}\\ & & +\displaystyle \frac{2i{\omega }_{D}^{s}{\mathcal{A}}{\mathcal{B}}}{\mu (\mu +1)}\left[{(\mu +1)}^{2}-\mu (\mu -1){\mathcal{B}}\right].\end{eqnarray} \tag{ 80 }$$

In the limit of vanishing magnetic fields (${{\mathcal{M}}}^{2}\to \mu$ and ${\mathcal{A}}\to {(1+\mu )}^{-1}$), ${{\mathcal{E}}}_{S}={{\mathcal{E}}}_{{\rm{M}}}=0$, and Equation (75) reduces to

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{s}}}\left({{\rm{\Xi }}}_{1}\right)}{{{ik}}_{y}} & = & \displaystyle \frac{{(\mu -1)}^{2}}{{\mu }^{2}}\displaystyle \frac{1}{{{\rm{\Sigma }}}_{0}^{s-}}\left({U}_{1}^{s-}+\displaystyle \frac{{{du}}_{T0}^{s-}}{{dx}}{Z}_{1}+i{\omega }_{D}^{s}{Z}_{1}\right)\\ & & -\displaystyle \frac{{q}_{c}{\rm{\Omega }}L}{{u}_{T0}{{\rm{\Sigma }}}_{0}^{s-}}{V}_{1}^{s-},\end{eqnarray} \tag{ 81 }$$

identical to Equation (A8) of Paper I. Paper I showed that the ${\omega }_{D}$^sZ₁ term in Equation (81) is responsible for the production of the PV at a disturbed shock, while the terms involving ${U}_{1}^{s-}$ and ${V}_{1}^{s-}$ suppress the PV due to shock compression and background shear, respectively.