text:
# Perturbation configurations in a two-fluid system of singular isothermal disks

## 1 Introduction

The theoretical model problem here is one that involves large-scale stationary density waves in a composite system of gravitationally coupled stellar and gaseous disks with the two fluid disks being both idealized as razor-thin singular isothermal disks (SIDs). For the gravitational effect of a background axisymmetric dark-matter halo, we also consider a background composite system of two coupled partial SIDs (Syer & Tremaine 1996; Shu, Laughlin, Lizano & Galli 2000; Lou 2002). We search for stationary coplanar perturbation configurations in such a composite SID system.

Over nearly four decades, there have been numerous theoretical studies on perturbation and stability properties of a composite system of stellar and gaseous disks, mostly in galactic contexts. The local dispersion relation for galactic spiral density waves in a composite disk system of stars and gas was first derived by Lin & Shu (1966, 1968), where the collective behavior of stars was described by a stellar distribution function while the gaseous disk was treated as an isothermal fluid. Their pioneer work was later followed up by Kato (1972), who examined oscillations and overstability of density waves in a similar formalism. Using the two-fluid approach, Jog & Solomon (1984a,b) studied the growth of local axisymmetric perturbations in gravitationally coupled stellar and gaseous disks. They found that a composite disk system can be unstable owing to the gravitational coupling, even though the stellar and gaseous disks can be separately stable. Bertin & Romeo (1988) studied the role of gas on global spiral modes in a two-fluid model framework. Vandervoort (1991a,b) investigated effects of interstellar gas on oscillations and stability of spheroidal galaxies. Romeo (1992) considered the stability of a two-component disk of finite thickness. The two-fluid formalism has been adopted into a modal analysis for morphologies of disk galaxies (Lowe et al. 1994). Elmegreen (1995) and Jog (1996) simultaneously approached a similar stability problem to derive an effective $`Q_{eff}`$ parameter (Safronov 1960; Toomre 1964) for axisymmetric two-fluid instabilities relevant to a disk galaxy. Lou & Fan (1998b) explored basic properties of open and tight-winding spiral density-wave modes in a composite disk system using the two-fluid formalism.

Since Mestel (1963), the concept of SIDs has attracted considerable theoretical interests in various contexts of disk dynamics in general (Zang 1976; Toomre 1977; Lemos, Kalnajs & Lynden-Bell 1991; Syer & Tremaine 1996; Lynden-Bell & Lemos 1999; Goodman & Evans 1999; Shu et al. 2000; Galli et al. 2001). These studies provide useful information for the research of star-formation (e.g. Shu et al. 1999 on formation and collapse of cloud cores in the birth of stars and planetary systems) and galactic structure communities (e.g., Bertin & Lin 1996 on the structure of barred and spiral galaxies and Crane et al. 1993 on the light cusps seen in the nuclei of galaxies), and deepen our understanding for the dynamics of self-gravitating configurations that have a power-law density distribution. In particular, Shu et al. (2000) derived perturbation solutions and performed a stability analysis on an isopedically magnetized SID with a flat rotation curve. The background equilibrium surface mass density distribution was assumed to bear a power-law radial variation. They found both stationary aligned and unaligned logarithmic configurations111Aligned disturbances involve distorted streamlines with the maximum and minimum radii at different radial locations lined up in the azimuth, while for unaligned or spiral disturbances, distorted streamlines with the maximum and minimum radii shifted in azimuth at different radial locations (Kalnajs 1973). and offer physical interpretations for the marginal instability curves in a single SID. In contrast to their work on a single SID system, we are mainly interested in the situation of a disk system composed of a stellar SID and a gaseous SID. As it is more realistic to consider large-scale dynamics of stellar and gaseous disks in a disk galaxy, the investigation on such a composite SID system can reveal more useful information. Motivated by such a prospect and combining our prior experience of treating composite disk system (Lou & Fan 1998b, 2000a, b), we search for both stationary aligned and unaligned or spiral logarithmic configurations in a composite SID system and discuss their stability properties (Safronov 1960; Toomre 1964; Goldreich & Lynden-Bell 1965; Ostriker & Peebles 1973; Elmegreen 1995; Jog 1996; Shu et al. 2000).

Large-scale dynamics of the gas disk is a necessary component in the overall density-wave scenario. Moreover, most of the observational diagnostics involve processes of gaseous interstellar medium (ISM) on various sub-scales. In more realistic terms, numerical simulation experiments are indispensable for studying linear and nonlinear dynamical processes. For numerical code development in particular, the analytical solutions here not only offer important physical insights but also serve as valuable tools of benchmarking. Complementarily, numerical simulations starting from or based upon these stationary perturbation solutions can lead to insights for equilibrium, stability and nonlinear processes (e.g., spiral shocks to trigger star formation activities). Numerical simulations for clusters of galaxies and numerical simulations for galactic disk dynamics, although on totally different scales, are similar in several fundamental aspects, that is, both involve massive dark-matter halo, N-body gravitational interaction (galaxies as mass points for a cluster versus stars as mass points for a disk galaxy), and gas dynamics. The main difference lies in the geometry involved, namely, grossly spherical geometry for a typical galaxy cluster and a disk geometry for a spiral galaxy. Therefore, with proper adaptation, numerical codes designed for a galaxy cluster can be applied to a disk galaxy or vice versa.

For spiral galaxies, a more realistic model would involve a magnetized gas disk gravitationally coupled to a stellar disk with differential rotation in the presence of a massive dark-matter halo. For a single magnetized SID (MSID) with a coplanar magnetic field, we have recently derived solutions for stationary aligned and logarithmic spiral configurations (Lou 2002) and pointed out slow MSID configurations can persist in an extended radial range of a disk with flat rotation curve (Lou & Fan 2002) as in the case of NGC 6946 for a spiral pattern of interlaced optical and magnetic field arms (e.g., Fan & Lou 1996; Lou & Fan 1998; Frick et al. 2000). The present hydrodynamic problem is not only interesting by itself, but also serves as a necessary step for constructing more realistic stationary configurations in a composite system of a stellar SID and an MSID.

The paper is structured as follows. In Section 2, we present the basic two-fluid formalism for the composite system of SIDs coupled by self-gravity, derive equilibrium properties for both SIDs, and obtain the linearized equations for small disturbances. Aligned and unaligned solutions for stationary perturbations in a composite SID system and their stability properties are studied and discussed in Section 3. We explore various parameter regimes for physical solutions and compare them with the results of the single SID case. The analysis is extended to a composite partial SID system in Section 4. We summarize computational results and analysis in Section 5 and discuss potential galactic applications. Specific details are included in Appendices A$``$D.

## 2 Two-fluid SIDs formalism

For a model containing sufficient physics and for mathematical simplicity, we adopt the two-fluid formalism for large-scale stationary aligned or unaligned disturbances in a background rotational equilibrium with axisymmetry. In dealing with singularities and resonances, it would be physically more accurate to adopt the formalism of distribution functions (e.g. Lin & Shu 1966; Julian & Toomre 1966; Binney & Tremaine 1987). For the purpose of this study, such irregularities and resonances do not arise and the two-fluid approach will offer useful information for us to learn. In this section, we provide the basic equations for the two-fluid system, composed of a stellar disk and a gaseous disk. Given qualifications and assumptions, equilibrium properties of the stellar and gaseous SIDs with flat rotation curves (allowed to be different in a consistent manner) are summarized. We derive coplanar perturbation equations in both stellar and gaseous SIDs, respectively.

### 2.1 Two sets of coupled fluid equations

For expediency, the two SIDs located at $`z=0`$ are treated as infinitesimally thin, which are sometimes referred to as razor-thin disks. For physical variables, we shall use superscripts and/or subscripts $`s`$ to denote the stellar disk and superscripts and/or subscripts $`g`$ to denote the gaseous disk. The two razor-thin rotating disks in a composite system are modelled as two fluids coupled through the mutual gravitational interaction. In the present formulation of large-scale perturbations, we ignore nonideal effects such as viscosity, resistivity, and thermal diffusion etc. Then the two fully nonlinear equation sets for the stellar and gaseous disks can be written out using cylindrical coordinates (r,$`\phi ,z)`$ in the $`z=0`$ plane. For the stellar disk, we have

$$\frac{\mathrm{\Sigma }^s}{t}+\frac{1}{r}\frac{}{r}(r\mathrm{\Sigma }^su^s)+\frac{1}{r^2}\frac{}{\phi }(\mathrm{\Sigma }^sj^s)=0,$$
(1)
$$\frac{u^s}{t}+u^s\frac{u^s}{r}+\frac{j^s}{r^2}\frac{u^s}{\phi }\frac{j^{s2}}{r^3}=\frac{1}{\mathrm{\Sigma }^s}\frac{}{r}(a_s^2\mathrm{\Sigma }^s)\frac{\varphi }{r},$$
(2)
$$\frac{j^s}{t}+u^s\frac{j^s}{r}+\frac{j^s}{r^2}\frac{j^s}{\phi }=\frac{1}{\mathrm{\Sigma }^s}\frac{}{\phi }(a_s^2\mathrm{\Sigma }^s)\frac{\varphi }{\phi }.$$
(3)
In parallel, we have for the gaseous disk

$$\frac{\mathrm{\Sigma }^g}{t}+\frac{1}{r}\frac{}{r}(r\mathrm{\Sigma }^gu^g)+\frac{1}{r^2}\frac{}{\phi }(\mathrm{\Sigma }^gj^g)=0,$$
(4)
$$\frac{u^g}{t}+u^g\frac{u^g}{r}+\frac{j^g}{r^2}\frac{u^g}{\phi }\frac{j^{g2}}{r^3}=\frac{1}{\mathrm{\Sigma }^g}\frac{}{r}(a_g^2\mathrm{\Sigma }^g)\frac{\varphi }{r},$$
(5)
$$\frac{j^g}{t}+u^g\frac{j^g}{r}+\frac{j^g}{r^2}\frac{j^g}{\phi }=\frac{1}{\mathrm{\Sigma }^g}\frac{}{\phi }(a_g^2\mathrm{\Sigma }^g)\frac{\varphi }{\phi }.$$
(6)
The coupling of the two sets of fluid equations is due to the gravitational potential through Poisson’s integral

$`\varphi (r,\phi ,t)={\displaystyle 𝑑\psi _0^{\mathrm{}}\frac{G\mathrm{\Sigma }(r^{},\psi ,t)r^{}dr^{}}{\left[r^2+r^22rr^{}\mathrm{cos}(\psi \phi )\right]^{1/2}}},`$ (7)
where $`\mathrm{\Sigma }=\mathrm{\Sigma }^s+\mathrm{\Sigma }^g`$ is the total surface mass density. In equations $`(1)(7)`$, $`\mathrm{\Sigma }^s`$ is the stellar surface mass density, $`u^s`$ is the radial component of the fluid velocity, $`j^s`$ is the $`z`$component of the specific angular momentum, and $`a_s`$ is the stellar velocity dispersion (or an effective “isothermal sound speed”), $`a_s^2\mathrm{\Sigma }^s`$ stands for an effective pressure in the polytropic approximation, $`\varphi `$ is the total gravitational potential expressed in terms of Poisson’s integral. For physical variables of the gaseous disk, we simply replace the superscript $`s`$ by $`g`$ systematically. Here, we assume that the stellar and gaseous disks interact mainly through the mutual gravitational coupling on large scales (Jog & Solomon 1984a,b; Bertin & Romeo 1988; Romeo 1992; Elmegreen 1995; Jog 1996; Lou & Fan 1998b, 2000a,b).

### 2.2 Properties of an axisymmetric equilibrium

Before a coplanar perturbation analysis, one needs to adopt a background rotational equilibrium for the composite SID system consistent with equations $`(1)(7)`$. We assume axisymmetric background SIDs for both stellar and gaseous disks, with the same form of power-law surface mass densities ($`\mathrm{\Sigma }r^1`$) yet with different flat rotation curves. The equilibrium properties of the composite system can then be derived from the basic equations of Section 2.1. In applications, the divergence towards $`r0`$ may be bypassed by introducing a gradual transition from disk to bulge or artificial inner cut-outs. In theoretical analyses, such a $`r0`$ divergence poses a challenge of understanding the SID stability properties (Zang 1976; Toomre 1977; Lynden-Bell & Lemos 1999; Evans & Read 1998; Goodman & Evans 1999; Shu et al. 2000). In particular, Shu et al. (2000) have made a systematic investigation on the problem in an attempt to summarize and clarify the relevant issues and to resolve the discrepancies between the disparate viewpoints.

Using equations (2) and (5) for the background equilibrium with $`u_0^s=u_0^g=0`$, $`\mathrm{\Omega }_s=j_0^s/r^2`$, and $`\mathrm{\Omega }_g=j_0^g/r^2`$, we readily obtain

$$\begin{array}{c}\hfill \mathrm{\Sigma }_0^s=a_s^2(1+D_s^2)/[2\pi Gr(1+\delta )],\\ \hfill \mathrm{\Sigma }_0^g=a_g^2(1+D_g^2)\delta /[2\pi Gr(1+\delta )],\end{array}$$
(8)
where $`\delta \mathrm{\Sigma }_0^g/\mathrm{\Sigma }_0^s`$ is the background surface mass density ratio. Furthermore, one can write

$$\begin{array}{c}\hfill \mathrm{\Omega }_s=a_sD_s/r,\\ \hfill \mathrm{\Omega }_g=a_gD_g/r,\end{array}$$
(9)
$$\begin{array}{c}\hfill \kappa _s\{(2\mathrm{\Omega }_s/r)[d(r^2\mathrm{\Omega }_s)/dr]\}^{1/2}=\sqrt{2}\mathrm{\Omega }_s,\\ \hfill \kappa _g\{(2\mathrm{\Omega }_g/r)[d(r^2\mathrm{\Omega }_g)/dr]\}^{1/2}=\sqrt{2}\mathrm{\Omega }_g,\end{array}$$
(10)
$$a_s^2(D_s^2+1)=a_g^2(D_g^2+1),$$
(11)
where $`\mathrm{\Omega }_s`$ and $`\mathrm{\Omega }_g`$ are the mean angular rotation speeds of the stellar and gaseous SIDs, respectively, $`\kappa _s`$ and $`\kappa _g`$ are the corresponding epicyclic frequencies, and $`D_s`$ and $`D_g`$ are the dimensionless parameters characterizing the level of rotation for the stellar and gaseous disks, respectively. The surface mass densities $`\mathrm{\Sigma }_s`$ and $`\mathrm{\Sigma }_g`$ of both stellar and gaseous SIDs take the power-law form of $`r^1`$. Note that condition (11) (also valid for a composite system of two coupled partial SIDs discussed later) is derived from the equilibrium radial momentum equations (2) and (5) using the polytropic approximation for both SIDs. It implies an intimate relation among the two SID rotation speeds, the stellar velocity dispersion and the gas sound speed. In this aspect, it is different from the usual prescription for a composite disk system (e.g., Jog & Solomon 1984a, b; Bertin & Romeo 1988; Elmegreen 1995; Jog 1996; Lou & Fan 1998b). The two SID rotation speeds will be different as long as $`a_sa_g`$. This may induce streaming instabilities when the difference between the two SID rotation speeds is sufficiently large.

Here, we introduce two useful parameters for a composite system of two gravitationally coupled SIDs. The first one is the SID surface mass density ratio $`\delta \mathrm{\Sigma }_0^g/\mathrm{\Sigma }_0^s`$. For late-type mature spiral galaxies, gas materials are less than the stellar mass, that is, $`\delta <1`$. However, for primordial disk galaxies, the gas materials exceed the stellar mass in general. Therefore, in our analysis and computations, both cases of $`\delta <1`$ and $`\delta 1`$ are considered. The second parameter is $`\beta a_s^2/a_g^2`$ for the square of the ratio of the velocity dispersion of stellar disk to the sound speed of gas disk. Typically, the stellar velocity dispersion is larger than the sound speed of gas disk. We therefore take $`a_s^2>a_g^2`$ or $`\beta >1`$. By “isothermal”, we mean constant $`a_s^2`$ and constant $`a_g^2`$. For actual spiral galaxies, this represents a gross simplification. In the special situation of $`a_s^2=a_g^2`$, it follows from condition (11) that $`D_s^2=D_g^2`$. Thus, the two SIDs may be treated as a single SID, because the stellar and gas disks rotate with the same speed. We will show later that this $`\beta =1`$ case gives some familiar results of a single SID which can be regarded as the limiting case of a two-SID system. In our analysis, we shall only consider the case222Background equilibrium condition (11) guarantees that $`D_g^20`$ as long as $`D_s^20`$ when $`\beta 1`$. of $`\beta 1`$.

### 2.3 Coplanar perturbation equations in SIDs

We now generally consider small nonaxisymmetric perturbations, marked along a physical variable with a subscript 1, in both stellar and gaseous SIDs. For example,

$$\mathrm{\Sigma }^s=\mathrm{\Sigma }{}_{0}{}^{s}+\mathrm{\Sigma }{}_{1}{}^{s},\mathrm{\Sigma }^g=\mathrm{\Sigma }{}_{0}{}^{g}+\mathrm{\Sigma }{}_{1}{}^{g},$$
(12)
$$u^s=u_0^s+u_1^s,u^g=u_0^g+u_1^g,$$
(13)
$$j^s=j_0^s+j_1^s,j^g=j_0^g+j_1^g,$$
(14)
$$\mathrm{\Sigma }=\mathrm{\Sigma }^s+\mathrm{\Sigma }^g=(\mathrm{\Sigma }{}_{0}{}^{s}+\mathrm{\Sigma }{}_{0}{}^{g})+(\mathrm{\Sigma }{}_{1}{}^{s}+\mathrm{\Sigma }{}_{1}{}^{g}).$$
(15)
Substituting expressions $`(12)(15)`$ into full equations $`(1)(7)`$ and linearizing about the axisymmetric background denoted by a subscript 0, we derive

$$\frac{\mathrm{\Sigma }_1^s}{t}+\frac{1}{r}\frac{}{r}(r\mathrm{\Sigma }_0^su_1^s)+\mathrm{\Omega }_s\frac{\mathrm{\Sigma }_1^s}{\phi }+\frac{\mathrm{\Sigma }_0^s}{r^2}\frac{j_1^s}{\phi }=0,$$
(16)
$$\frac{u_1^s}{t}+\mathrm{\Omega }_s\frac{u_1^s}{\phi }2\mathrm{\Omega }_s\frac{j_1^s}{r}=\frac{}{r}\left(a_s^2\frac{\mathrm{\Sigma }_1^s}{\mathrm{\Sigma }_0^s}+\varphi _1\right),$$
(17)
$$\frac{j_1^s}{t}+r\frac{\kappa _s^2}{2\mathrm{\Omega }_s}u_1^s+\mathrm{\Omega }_s\frac{j_1^s}{\phi }=\frac{}{\phi }\left(a_s^2\frac{\mathrm{\Sigma }_1^s}{\mathrm{\Sigma }_0^s}+\varphi _1\right)$$
(18)
for the stellar SID, and

$$\frac{\mathrm{\Sigma }_1^g}{t}+\frac{1}{r}\frac{}{r}(r\mathrm{\Sigma }_0^gu_1^g)+\mathrm{\Omega }_g\frac{\mathrm{\Sigma }_1^g}{\phi }+\frac{\mathrm{\Sigma }_0^g}{r^2}\frac{j_1^g}{\phi }=0,$$
(19)
$$\frac{u_1^g}{t}+\mathrm{\Omega }_g\frac{u_1^g}{\phi }2\mathrm{\Omega }_g\frac{j_1^g}{r}=\frac{}{r}\left(a_g^2\frac{\mathrm{\Sigma }_1^g}{\mathrm{\Sigma }_0^g}+\varphi _1\right),$$
(20)
$$\frac{j_1^g}{t}+r\frac{\kappa _g^2}{2\mathrm{\Omega }_g}u_1^g+\mathrm{\Omega }_g\frac{j_1^g}{\phi }=\frac{}{\phi }\left(a_g^2\frac{\mathrm{\Sigma }_1^g}{\mathrm{\Sigma }_0^g}+\varphi _1\right)$$
(21)
for the gaseous SID, with the total gravitational potential perturbation given by

$`\varphi _1(r,\phi ,t)={\displaystyle 𝑑\psi _0^{\mathrm{}}\frac{G(\mathrm{\Sigma }_1^s+\mathrm{\Sigma }_1^g)r^{}dr^{}}{\left[r^2+r^22rr^{}\mathrm{cos}(\psi \phi )\right]^{1/2}}}.`$ (22)
Assuming a Fourier periodic form of $`\mathrm{exp}[i(\omega tm\phi )]`$ for perturbation solutions in general (after taking the real part), we write for coplanar perturbations in the stellar disk as

$$\begin{array}{c}\hfill \mathrm{\Sigma }{}_{1}{}^{s}=\mu ^s(r)\mathrm{exp}[i(\omega tm\phi )],\\ \hfill u_1^s=U^s(r)\mathrm{exp}[i(\omega tm\phi )],\\ \hfill j_1^s=J^s(r)\mathrm{exp}[i(\omega tm\phi )],\end{array}$$
(23)
for coplanar perturbations in the gaseous disk as

$$\begin{array}{c}\hfill \mathrm{\Sigma }{}_{1}{}^{g}=\mu ^g(r)\mathrm{exp}[i(\omega tm\phi )],\\ \hfill u_1^g=U^g(r)\mathrm{exp}[i(\omega tm\phi )],\\ \hfill j_1^g=J^g(r)\mathrm{exp}[i(\omega tm\phi )],\end{array}$$
(24)
and for the total gravitational potential perturbation as

$$\varphi _1=V(r)\mathrm{exp}[i(\omega tm\phi )]$$
(25)
in the SID plane at $`z=0`$. By substituting expressions $`(23)(25)`$ into equations $`(16)(22)`$, we derive for the stellar disk

$$i(\omega m\mathrm{\Omega }_s)\mu ^s+\frac{1}{r}\frac{d}{dr}(r\mathrm{\Sigma }_0^sU^s)im\mathrm{\Sigma }{}_{0}{}^{s}\frac{J^s}{r^2}=0,$$
(26)
$$i(\omega m\mathrm{\Omega }_s)U^s2\mathrm{\Omega }_s\frac{J^s}{r}=\frac{d}{dr}\left(a_s^2\frac{\mu ^s}{\mathrm{\Sigma }_0^s}+V\right),$$
(27)
$$i(\omega m\mathrm{\Omega }_s)J^s+r\frac{\kappa _s^2}{2\mathrm{\Omega }_s}U^s=im\left(a_s^2\frac{\mu ^s}{\mathrm{\Sigma }_0^s}+V\right),$$
(28)
for the gaseous disk

$$i(\omega m\mathrm{\Omega }_g)\mu ^g+\frac{1}{r}\frac{d}{dr}(r\mathrm{\Sigma }_0^gU^g)im\mathrm{\Sigma }{}_{0}{}^{g}\frac{J^g}{r^2}=0,$$
(29)
$$i(\omega m\mathrm{\Omega }_g)U^g2\mathrm{\Omega }_g\frac{J^g}{r}=\frac{d}{dr}\left(a_g^2\frac{\mu ^g}{\mathrm{\Sigma }_0^g}+V\right),$$
(30)
$$i(\omega m\mathrm{\Omega }_g)J^g+r\frac{\kappa _g^2}{2\mathrm{\Omega }_g}U^g=im\left(a_g^2\frac{\mu ^g}{\mathrm{\Sigma }_0^g}+V\right),$$
(31)
and for the total gravitational potential perturbation

$$V(r)=𝑑\psi _0^{\mathrm{}}\frac{G(\mu ^s+\mu ^g)r^{}dr^{}}{\left[r^2+r^22rr^{}\mathrm{cos}\psi \right]^{1/2}}.$$
(32)
We now use equations (27) and (28) to express $`U^s`$ and $`J^s`$ in terms of $`\mathrm{\Psi }^sa_s^2\mu ^s/\mathrm{\Sigma }_0^s+V`$ for the stellar SID and similarly, use equations (30) and (31) to express $`U^g`$ and $`J^g`$ in terms of $`\mathrm{\Psi }^ga_g^2\mu ^g/\mathrm{\Sigma }_0^g+V`$ for the gaseous SID. The resulting expressions then become

$$U^s=\frac{i}{(\omega m\mathrm{\Omega }_s)^2\kappa _s^2}\left[2\mathrm{\Omega }_s\frac{m}{r}+(\omega m\mathrm{\Omega }_s)\frac{d}{dr}\right]\mathrm{\Psi }^s$$
(33)
and

$$\frac{J^s}{r}=\frac{1}{(\omega m\mathrm{\Omega }_s)^2\kappa _s^2}\left[(\omega m\mathrm{\Omega }_s)\frac{m}{r}\frac{\kappa _s^2}{2\mathrm{\Omega }_s}\frac{d}{dr}\right]\mathrm{\Psi }^s$$
(34)
for the stellar SID, and

$$U^g=\frac{i}{(\omega m\mathrm{\Omega }_g)^2\kappa _g^2}\left[2\mathrm{\Omega }_g\frac{m}{r}+(\omega m\mathrm{\Omega }_g)\frac{d}{dr}\right]\mathrm{\Psi }^g$$
(35)
and

$$\frac{J^g}{r}=\frac{1}{(\omega m\mathrm{\Omega }_g)^2\kappa _g^2}\left[(\omega m\mathrm{\Omega }_g)\frac{m}{r}\frac{\kappa _g^2}{2\mathrm{\Omega }_g}\frac{d}{dr}\right]\mathrm{\Psi }^g$$
(36)
for the gaseous SID, respectively.

Substitution of expressions (33) and (34) into equation (26) leads to

$`0=(\omega m\mathrm{\Omega }_s)\mu ^s+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}`$
$`\times \left\{{\displaystyle \frac{r\mathrm{\Sigma }_0^s}{(\omega m\mathrm{\Omega }_s)^2\kappa _s^2}}\left[2\mathrm{\Omega }_s{\displaystyle \frac{m}{r}}+(\omega m\mathrm{\Omega }_s){\displaystyle \frac{d}{dr}}\right]\mathrm{\Psi }^s\right\}`$
$`{\displaystyle \frac{m\mathrm{\Sigma }_0^s}{r[(\omega m\mathrm{\Omega }_s)^2\kappa _s^2]}}\left[(\omega m\mathrm{\Omega }_s){\displaystyle \frac{m}{r}}{\displaystyle \frac{\kappa _s^2}{2\mathrm{\Omega }_s}}{\displaystyle \frac{d}{dr}}\right]\mathrm{\Psi }^s`$ (37)
for the stellar SID. Similarly, substitution of expressions (35) and (36) into equation (29) leads to

$`0=(\omega m\mathrm{\Omega }_g)\mu ^g+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}`$
$`\times \left\{{\displaystyle \frac{r\mathrm{\Sigma }_0^g}{(\omega m\mathrm{\Omega }_g)^2\kappa _g^2}}\left[2\mathrm{\Omega }_g{\displaystyle \frac{m}{r}}+(\omega m\mathrm{\Omega }_g){\displaystyle \frac{d}{dr}}\right]\mathrm{\Psi }^g\right\}`$
$`{\displaystyle \frac{m\mathrm{\Sigma }_0^g}{r[(\omega m\mathrm{\Omega }_g)^2\kappa _g^2]}}\left[(\omega m\mathrm{\Omega }_g){\displaystyle \frac{m}{r}}{\displaystyle \frac{\kappa _g^2}{2\mathrm{\Omega }_g}}{\displaystyle \frac{d}{dr}}\right]\mathrm{\Psi }^g`$ (38)
for the gaseous SID. Equations (37) and (38) are to be solved with Poisson’s integral (32).

For stationary solutions ($`\omega =0`$) with zero pattern speed, equations (37) and (38) above can be cast into the forms of equations (39) and (40) below by invoking the background equilibrium conditions $`(8)(11)`$. That is,

$`m[\mu ^s+{\displaystyle \frac{1}{D_s^2(m^22)}}({\displaystyle \frac{m^2}{r}}2{\displaystyle \frac{d}{dr}}r{\displaystyle \frac{d^2}{dr^2}})`$
$`\times (r\mu ^s+{\displaystyle \frac{1+D_s^2}{2\pi G}}{\displaystyle \frac{V}{1+\delta }})]=0`$ (39)
for the stellar SID, and

$`m[\mu ^g+{\displaystyle \frac{1}{D_g^2(m^22)}}({\displaystyle \frac{m^2}{r}}2{\displaystyle \frac{d}{dr}}r{\displaystyle \frac{d^2}{dr^2}})`$
$`\times (r\mu ^g+{\displaystyle \frac{1+D_g^2}{2\pi G}}{\displaystyle \frac{V\delta }{1+\delta }})]=0`$ (40)
for the gaseous SID. Again, equations (39) and (40) are to be solved simultaneously together with Poisson’s integral (32).

## 3 Aligned and unaligned cases

### 3.1 Aligned perturbation configurations

Let us now obtain the aligned stationary density wave patterns from equations (32), (39) and (40). We note in particular that aligned perturbations relate to purely azimuthal propagations of density waves (see Section 3.2 of Lou 2002).

#### 3.1.1 The $`|m|=0`$ Case: axisymmetric disturbances

We first examine the $`m=0`$ case. With $`\omega =m=0`$, the solution to equations $`(26)(32)`$ takes the forms of $`U^s=U^g=0`$, $`\mu ^s=K_1^s/r`$, $`\mu ^g=K_1^g/r`$, $`J^s=K_2^sr`$, $`J^g=K_2^gr`$, $`V=K\mathrm{ln}r`$, where the ratios of constants $`K_1^s/K`$, $`K_1^g/K`$, $`K_2^s/K`$, and $`K_2^g/K`$ are chosen such that equations (27), (30), (32) can be satisfied. However, such a “solution” merely represents a rescaling of one axisymmetric equilibrium to a neighboring axisymmetric equilibrium (Shu et al. 2000). This rescaling is allowed by equations (32), (39) and (40) but uninteresting in the present context. We now turn to cases with $`|m|1`$.

#### 3.1.2 Cases with $`|m|1`$: nonaxisymmetric disturbances

In power-law disks, we consider aligned perturbations that carry the same power-law dependence as the equilibrium SID does. By this assumption, we mean the following exact potential-density relation

$$\begin{array}{c}\hfill \mu _s=\sigma _s/r,\\ \hfill \mu _g=\sigma _g/r,\\ \hfill V=2\pi Gr\mu _s/|m|2\pi Gr\mu _g/|m|,\end{array}$$
(41)
where $`\sigma _s`$ and $`\sigma _g`$ are constant coefficients. It can be verified that Poisson’s integral (32) is satisfied (Shu et al. 2000). A direct substitution of equation (41) into equations (39) and (40) then leads to the following equations:

$$\mu ^s=\left(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2}\right)(H_1r\mu ^s+G_1r\mu ^g),$$
(42)
$$\mu ^g=\left(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2}\right)(H_2r\mu ^g+G_2r\mu ^s),$$
(43)
where coefficients $`H_1`$, $`H_2`$, $`G_1`$ and $`G_2`$ are defined by

$$\begin{array}{c}\hfill H_1\frac{1}{D_s^2(m^22)}\left(1\frac{D_s^2+1}{|m|}\frac{1}{1+\delta }\right),\\ \hfill H_2\frac{1}{D_g^2(m^22)}\left(1\frac{D_g^2+1}{|m|}\frac{\delta }{1+\delta }\right),\end{array}$$
(44)
$$\begin{array}{c}\hfill G_1\frac{D_s^2+1}{D_s^2(m^22)|m|}\frac{1}{1+\delta },\\ \hfill G_2\frac{D_g^2+1}{D_g^2(m^22)|m|}\frac{\delta }{1+\delta }..\end{array}$$
(45)
By substituting the forms of $`\mu ^s`$ and $`\mu ^g`$ given by equation (41) into equations (42) and (43) above, we readily obtain

$$(1H_1m^2)\mu ^s=G_1m^2\mu ^g$$
(46)
and

$$(1H_2m^2)\mu ^g=G_2m^2\mu ^s.$$
(47)
It then follows from equations (46) and (47) that

$$(1H_1m^2)(1H_2m^2)=G_1G_2m^4.$$
(48)
In the above procedure, we have assumed $`\delta 0`$ in order to eliminate $`\mu ^s`$ and $`\mu ^g`$ from both sides of equation (48). Otherwise the problem would reduce to that of a single SID.

Using definitions (44) and (45) for the expressions of $`H_1`$, $`H_2`$, $`G_1`$ and $`G_2`$, we obtain explicitly

$$\begin{array}{c}\hfill \left[1\frac{m^2}{D_s^2(m^22)}\left(1\frac{D_s^2+1}{|m|}\frac{1}{1+\delta }\right)\right]\\ \hfill \times \left[1\frac{m^2}{D_g^2(m^22)}\left(1\frac{D_g^2+1}{|m|}\frac{\delta }{1+\delta }\right)\right]\\ \hfill =\frac{D_s^2+1}{D_s^2(m^22)}\frac{D_g^2+1}{D_g^2(m^22)}\frac{m^2\delta }{(1+\delta )^2}\end{array}$$
(49)
which can be further simplified to

$$\begin{array}{c}\hfill \left[D_s^2(m^22)m^2\left(1\frac{D_s^2+1}{|m|}\frac{1}{1+\delta }\right)\right]\\ \hfill \times \left[D_g^2(m^22)m^2\left(1\frac{D_g^2+1}{|m|}\frac{\delta }{1+\delta }\right)\right]\\ \hfill =(D_s^2+1)(D_g^2+1)\frac{m^2\delta }{(1+\delta )^2}.\end{array}$$
(50)
By condition (11) for the background equilibrium, we have $`D_g^2=\beta (D_s^2+1)1`$ with $`\beta a_s^2/a_g^2`$. Equation (50) can be cast into a quadratic equation in terms of $`yD_s^2`$, namely

$$C_2y^2+C_1y+C_0=0,$$
(51)
where

$$C_2=\beta (m^22)(m^2+|m|2),$$
$`(51a)`$
$$\begin{array}{c}\hfill C_1=\left(2m^2+4\frac{2m\delta }{1+\delta }\frac{4m}{1+\delta }+\frac{2m^3}{1+\delta }\right)\beta \\ \hfill \text{ }+\frac{2m}{1+\delta }2m^4+6m^2\frac{2m^3}{1+\delta }4,(51b)\end{array}$$

$$\begin{array}{c}\hfill C_0=\left(\frac{m^3}{1+\delta }+2m^2m^4\frac{2m}{1+\delta }\frac{m^3\delta }{1+\delta }\right)\beta \\ \hfill +\frac{2m}{1+\delta }2m^2\frac{2m^3}{1+\delta }+2m^4.(51c)\end{array}$$

The physical meaning is that for aligned nonaxisymmetric stationary configurations to exist in a composite SID system, condition (50), or equivalently, condition (51) must be satisfied for proper values of $`D_s^2`$ given specified parameters $`m`$, $`\delta `$ and $`\beta `$. In essence, a purely azimuthal propagation of density wave in opposite direction relative to the SID rotation needs to be counterbalanced by disk rotation in order to appear stationary in an inertial frame of reference. There are two possible sets of density waves in a composite disk system (e.g., Lou & Fan 1998b). Therefore, there are two possible disk rotation speeds that may produce stationary perturbation configurations in general.

Mathematically, one can show that equation (51) always has two real solutions for $`yD_s^2`$, except for the special case of $`m=1`$. The two real solutions are not necessarily both physically valid due to the nonnegative requirement of $`D_s^20`$. Nevertheless, it can be proven later that there always exists at least one physical solution with $`yD_s^20`$.

For the aligned case of $`m=1`$, it turns out that $`C_2=C_1=C_0=0`$ by definitions $`(51a)(51c)`$ and $`\mu ^g/\mu ^s=\mathrm{\Sigma }_0^g/\mathrm{\Sigma }_0^s`$ according to equation (46). Equation (51) can therefore be satisfied for arbitrary values of $`D_s^2`$. This is quite similar to the $`m=1`$ case of a single SID studied by Shu et al. (2000). While noting earlier concerns of spurious artifacts of the analytical method about this result (Zang 1976; Toomre 1977), Shu et al (2000) expressed a strong belief that such aligned eccentric displacements are possible alternative states of equilibria for extended SIDs. Whether such $`m=1`$ case represents a trivial translation of the origin of coordinates, we note for example that contours of surface mass density (background plus perturbation ones) are given by $`\mathrm{\Sigma }=(\mathrm{\Sigma }_0r)/r+\sigma \mathrm{cos}\phi /r=`$constant according to equation (41), while contours of surface mass density resulting from a small shift $`a`$ along the $`x`$axis of the origin from $`x=0`$, $`y=0`$ to $`x=a`$, $`y=0`$ are given by $`\mathrm{\Sigma }=[(\mathrm{\Sigma }_0r)/r](1+a\mathrm{cos}\phi /r)=`$constant. It is then clear that $`m=1`$ aligned perturbations are not trivial translation of the origin of coordinates for the present model problem.333In analyzing stability properties of Maclaurin disks (Takahara 1976; Smith 1979), the marginal stationary case of $`m=1`$ represents a trivial translation of the origin of coordinates involving a uniform velocity perturbation as shown in footnote $`14`$ on p319 of Binney & Tremaine (1987). Using the same analytical criterion here, one can demonstrate that contours of the surface mass density (see eqns 5-99 and 5-117 of Binney & Tremaine 1987) are equivalent to a small shift of the origin of coordinates. On this ground, it seems plausible that stationary aligned eccentric $`|m|=1`$ configurations are possible alternative equilibria for a composite SID system, at least mathematically. We note, however, for a single partial SID (Lou 2002) and for a composite system of two coupled partial SIDs (discussed later), such stationary aligned eccentric $`|m|=1`$ configurations are not allowed. In other words, stationary aligned eccentric $`|m|=1`$ configurations might be possible for protostellar disks where dark-matter halos are not involved, but may not occur for disk galaxies where massive dark-matter halos are known to exist.

We now study cases of $`m2`$. As noted earlier, in a real spiral galaxy, the stellar velocity dispersion $`a_s`$ is usually larger than the sound speed $`a_g`$ of the gas disk444For example, one may take $`a_g=7\text{km s}^1`$ and $`a_s=30\text{km s}^1`$ for a typical late-type spiral galaxy., that is, $`\beta >1`$ in our model. Theoretically, the case of $`\beta =1`$ means that the stellar disk and gaseous disk have the same velocity dispersion, and thus the same rotation parameter $`D_s=D_g`$ by condition (11). In some sense, we may then treat the two SIDs as one single SID with $`\mathrm{\Sigma }=\mathrm{\Sigma }^s+\mathrm{\Sigma }^g`$. It is expected that the solution should have something in common with a single SID system (Shu et al. 2000). We consider below the case of $`\beta =1`$ for analytical solutions and for insights of general solution properties.

With $`\beta =1`$, equations (50) or (51) can be reduced to

$$(|m|1)[D_s^2(m^22)m^2][(|m|+2)D_s^2|m|]=0$$
(52)
for arbitrary $`\delta `$ values. Note that expressions $`(51a)`$, $`(51b)`$ and $`(51c)`$ are independent of $`\delta `$ when $`\beta =1`$. We therefore have two solutions $`D_s^2=m^2/(m^22)`$ and $`D_s^2=|m|/(|m|+2)`$ for $`|m|2`$, and arbitrary values of $`D_s^2`$ for $`|m|=1`$. Note that subsonic rotation $`D_s^2=|m|/(|m|+2)`$ for $`|m|2`$ is simply the result of a single SID (equation (27) of Shu et al. 2000), and the supersonic rotation $`D_s^2=m^2/(m^22)`$ for $`|m|2`$ is a novel feature of a composite SID system.

For the phase relationship between the two surface mass perturbations $`\mu ^s`$ and $`\mu ^g`$, we note that equation (46) can be written in the form of

$$\frac{\mu ^g}{\mu ^s}=\frac{1H_1m^2}{G_1m^2}=1\frac{[D_s^2(m^22)m^2](1+\delta )}{|m|(D_s^2+1)}.$$
(53)
For supersonic rotations with $`D_s^2=m^2/(m^22)`$, one simply has $`\mu ^g/\mu ^s=1`$, while for subsonic rotations with $`D_s^2=|m|/(|m|+2)`$, one has $`\mu ^g/\mu ^s=\delta `$ (Lou & Fan 1998b).

We have noted earlier (Lou & Fan 2002; Lou 2002) that stationary aligned nonaxisymmetric configurations in an inertial reference frame correspond to purely azimuthal propagations of density waves counterbalanced by SID rotation. Physically, $`\mu ^g=\mu ^s`$ simply means that surface mass density perturbations of stellar and gaseous SIDs are completely out of phase to reduce the effect of gravity. This in turn implies a faster azimuthal density wave speed and thus requires a faster disk rotation (larger $`D_s^2`$) to maintain a stationary configuration. The case of $`\mu ^g=\delta \mu ^s`$ means that mass disturbance in gaseous disk is in phase with the mass disturbance in stellar disk, and their ratio is the same as the ratio of the background surface mass densities of the two SIDs. As the effect of self-gravity is enhanced, the azimuthal density waves speed is slower and thus a slower disk rotation (smaller $`D_s^2`$) is needed to sustain a stationary configuration.

For more realistic situation of $`\beta >1`$ in general, we readily obtain two branches of solution for $`yD_s^2`$ from quadratic equation (51) as functions of $`\beta `$ when $`|m|2`$ and $`\delta `$ are specified. It is found that both solution branches (upper branch $`y_1`$ and lower branch $`y_2`$) monotonically decrease for increasing values of $`\beta `$, and asymptotically approach different limits when $`\beta `$ goes to infinity (see three numerical examples shown in Fig. 1 and Appendix B for the variation trends).

It can be shown analytically that the larger solution branch $`y_1`$ (i.e., the higher branch) satisfies the following inequality

$$\frac{m^2}{m^22}\frac{2m(m+1)}{(m^22)(m+2)(1+\delta )}<y_1<\frac{m^2}{m^22},$$
(54)
where the lower bound on the left-hand side, which is always positive, is obtained by taking the limit of $`\beta \mathrm{}`$ and the upper bound on the right-hand side is determined by simply taking $`\beta =1`$. One can also show analytically that the smaller solution branch $`y_2`$ (i.e., the lower branch) satisfies another inequality

$$1<y_2<|m|/(|m|+2),$$
(55)
where the lower bound on the left-hand side is derived by letting $`\beta \mathrm{}`$ and the upper bound on the right is obtained by simply taking $`\beta =1`$. Note that the left-hand side of $`y_1`$ in inequality (54) is always greater than the right-hand side of $`y_2`$ in inequality (55), so that the two solution branches will not intersect with each other and $`y_1`$ remains always greater than $`y_2`$. In other words, for a given $`|m|`$, while there are two possible aligned solutions at maximum, there is only one aligned configuration that a composite SID system can support at a time. However, this does not exclude the possibility that for different values of $`|m|`$, more than one aligned configurations may be sustained in a composite SID system simultaneously.

By inequalities (54) and (55) and by numerical examples shown in Figure 1, the lower branches are subsonic while the upper branches are positive but may vary from supersonic ($`\beta 1`$) to subsonic for sufficiently large $`\beta `$ values when $`\delta `$ is small enough (see the case of $`m=2`$ and $`\delta =0.2`$). For late-type spiral galaxies, $`\delta `$ ranges from 0.05 to 0.1 and $`\beta `$ may take values between $`1520`$. Based on the variation trend displayed in Fig. 1, it is then possible for late-type disk galaxies to support stationary bar configurations of the upper branch with subsonic SID rotations. It is also interesting to infer that for early-type young disk galaxies with large values of $`\delta `$, stationary bar configurations of the upper branch may be sustained by supersonic SID rotations.

By examples of Fig. 1, it also becomes clear that for cases of $`|m|2`$, the lower solution branch $`y_2`$ of equation (51) may become negative when $`\beta `$ is larger than a critical value $`\beta _c`$ for given values of $`m`$ and $`\delta `$. In other words, there is only one solution branch for $`\beta >\beta _c`$, i.e. the upper $`y_1=D_s^2`$ branch, which is physically meaningful; being negative, the lower $`y_2=D_s^2`$ branch becomes unphysical and should be discarded. This critical value $`\beta _c`$ can be determined analytically in terms of $`\delta `$ and $`m`$ as

$$\beta _c=\frac{2(m+1)}{(m+2)}\left[1+\frac{m}{(m^2+2m)\delta +m^22}\right].$$
(56)
For a given value of $`m`$, the critical value of $`\beta _c`$ decreases with increasing values of $`\delta `$. When $`\delta `$ approaches infinity, the critical value $`\beta _c`$ goes to a limiting value

$$\beta _{cLim}=2(m+1)/(m+2).$$
(57)
In contexts of disk galaxies, this seems to suggest that aligned barred configurations with $`m=2`$ of the lower solution branch may be sustained by subsonic SID rotations only when $`\beta <3/2`$ which is rather restrictive.

We now examine the phase relationship between the surface mass density perturbations $`\mu ^g`$ and $`\mu ^s`$ in general. One can demonstrate that for the upper $`y_1=D_s^2`$ branch, the inequality

$$1<\frac{\mu ^g}{\mu ^s}<\frac{|m|\delta }{(m^2+|m|2)\delta +m^22}$$
(58)
holds. It is found that $`\mu ^g/\mu ^s`$ ratio decreases with increasing values of $`D_s^2`$ (see Appendix C) and therefore increases with increasing values of $`\beta `$. The lower bound on the left-hand side of inequality (58) is determined by taking $`\beta =1`$ and the upper bound on the right-hand side of inequality (58) is obtained by letting $`\beta \mathrm{}`$. For the lower $`y_2=D_s^2`$ branch, we have the following inequality

$$\delta <\mu ^g/\mu ^s<|m|(1+\delta )1$$
(59)
for perturbation surface mass density ratio $`\mu ^g/\mu ^s`$. As the mass ratio $`\mu ^g/\mu ^s`$ also increases with increasing values of $`\beta `$, the lower bound on the left-hand side of inequality (59) is determined by taking $`\beta =1`$. However, the upper bound on the right-hand side of inequality (59) is determined according to the possible maximum value of $`\beta `$ that makes $`D_s^2`$ physically meaningful, i.e. the critical value $`\beta _c`$ given by expression (56). It is now clear that the $`D_s^2=y_1`$ solution branch is always characterized by a negative $`\mu ^g/\mu ^s`$ and the $`D_s^2=y_2`$ solution branch is always characterized by a positive $`\mu ^g/\mu ^s`$. In this regard, the $`\beta =1`$ case which can be analyzed thoroughly shows a general solution property.

As observational diagnostics, our analysis suggests two different types of stationary aligned bar configurations in terms of surface mass density distributions. For stationary bar configurations of the upper branch solution, gas bars or young stellar bars should be out of phase relative to old stellar bars. For stationary bar configurations of the lower branch solution, gas bars or young stellar bars should overlap with old stellar bars. Another point of interest is that for late-type spiral galaxies, $`\beta `$ may be as large as $`1520`$. By Fig.1, the lower $`y_2`$ solution branches are excluded and the secular barlike instabilities more likely occur along the upper $`y_2`$ solution branches.

In general, the two ratios $`\beta `$ and $`\delta `$ are independent of each other. When $`\delta `$ and $`\beta `$ are specified, equation (51) can be solved for two values of $`D_s^2`$ with different values of $`|m|`$. We now take on a specific numerical example below (see Fig. 1). For $`|m|=2`$ and $`\delta =5`$, equation (51) becomes

$$8\beta y^2+(146\beta )y14\beta +22=0,$$
(60)
which has two solutions for $`D_s^2=y`$, namely,

$$y_1=[7+3\beta +(49134\beta +121\beta ^2)^{1/2}]/(8\beta )$$
(61)
for $`1<\mu ^g/\mu ^s<5/11`$, and

$$y_2=[7+3\beta (49134\beta +121\beta ^2)^{1/2}]/(8\beta )$$
(62)
for $`5<\mu ^g/\mu ^s<11`$. For any value $`\beta >1`$, $`y_1`$ given by equation (61) is always positive. For $`1\beta 11/7`$, $`y_2`$ given by equation (62) is nonnegative, while when $`\beta >11/7`$, $`y_2`$ becomes negative and thus unphysical. Besides other constraints, this may imply that in a disk galaxy where the stellar velocity dispersion is much greater than the sound speed of the gaseous disk, the $`\mu ^g/\mu ^s<0`$ branch may tend to manifest. Note that the range of $`\mu ^g/\mu ^s`$ is also fully determined by inequality (58) once $`|m|`$ and $`\delta `$ are known.

For an numerical example of a late-type disk galaxy, we follow the similar procedure shown above and take parameters $`m=2`$, $`\delta =0.05`$ and $`\beta =4`$. Now equation (51) yields two solutions, namely, $`y_1=0.6008`$ and $`y_2=0.2972`$. Apparently, the second $`y_2`$ solution should be discarded since $`\beta `$ exceeds the critical value $`\beta _c`$ given by equation (56), i.e., $`\beta _c=11/4`$. The $`y_1`$ solution is physically meaningful with $`\mu ^g/\mu ^s=0.0822`$. In this case, stationary aligned perturbations in gas surface mass density and stellar surface mass density are out of phase with each other.

To identify the nature of the aligned solution condition (50) or (51), we now examine the closely relevant case studied by Shu et al. (2000) and reemphasize the perspective that stationary aligned perturbation configurations should be regarded as purely azimuthal propagation of hydrodynamic density waves (Lou 2002; Lou & Fan 2002). For this purpose, we write solution condition (50) in a physically suggestive form of

$`[\mathrm{\Omega }_s^2(m^22)m^2a_s^2/r^2+2\pi G\mathrm{\Sigma }_s^0|m|/r]`$
$`\times [\mathrm{\Omega }_g^2(m^22)m^2a_g^2/r^2+2\pi G\mathrm{\Sigma }_g^0|m|/r]`$
$`=4\pi ^2G^2\mathrm{\Sigma }_0^s\mathrm{\Sigma }_0^gm^2/r^2.`$ (63)
The right-hand side of equation (63) represents the mutual gravitational coupling between the two SIDs. In the absence of this coupling, the two factors on the left-hand side would give rise two separate conditions for stationary aligned perturbation configurations with $`|m|2`$ for stellar and gaseous SIDs, respectively. For a single SID, be it stellar or gaseous, equation (63) therefore reduces to the form of

$$m^2\mathrm{\Omega }^2=\kappa ^2+m^2a^2/r^22\pi G\mathrm{\Sigma }_0|m|/r.$$
(64)
According to the well-known dispersion relation of density waves derived under the tight-winding or WKBJ approximation (Lin & Shu 1964, 1966 or equation (39) of Shu et al. 2000), equation (64) can be recovered by replacing the radial wavenumber $`|k|`$ with the azimuthal wavenumber $`|m|/r`$ and setting $`\omega =0`$ in an inertia frame of reference as noted by Lou (2002) in the study of stationary MHD perturbation configurations in a single MSID. By this procedure, it is quite clear that equation (64) describes an azimuthal propagation of hydrodynamic density waves, and a stationary pattern in the sidereal frame of reference requires specific values of $`D^2`$ for different $`|m|`$ values (Lou 2002; Lou & Fan 2002). In reference to dispersion relation (3.20) of Lou & Fan (1998b) with a proper adjustment of two different rotation speeds of SIDs, it is also transparent that equation (63) represents a purely azimuthal propagation of density waves in a composite SID system. The two coupled SIDs rotate in such a way to render the azimuthal density wave pattern stationary in a sidereal frame of reference.

In the context of a single SID, Shu et al. (2000) suggested that stationary condition (64) represents the onset of bifurcations of axisymmetric SIDs to nonaxisymmetrical SIDs that are more centrally condensed when the rotation rate $`D`$ is systematically increased. This is analogous to classical Maclaurin-Jacobi or Maclaurin-Dedekind bifurcations from spheroids to ellipsoids (Chandrasekhar 1969; Tassoul 1978) and involve secular instabilities that may be induced by viscous dissipations or by gravitational radiation (Bardeen et al. 1977). Shu et al. (2000) indicated that the instability mechanisms for aligned and spiral instabilities are fundamentally different and emphasized that wave propagation plays no role in aligned perturbations in contrast to spiral perturbations (see discussions in their subsection 3.1). Base on our analysis, we would like to clarify here that both types of perturbations involve propagations of density waves, even though a radial wave propagation for the spiral case may lead to different instability mechanisms than those without radial wave propagation for the aligned case. It is also important to note that for stationary configurations to appear in a sidereal reference frame, these density waves travel in opposite sense of SID rotation relative to the SID system and may strike a balance between wave propagation and SID rotation.

#### 3.1.3 Secular Barlike Instabilities

Strictly speaking, we have constructed analytically stationary aligned nonaxisymmetric configurations in a composite SID system. Are they stable or do they merely represent transition states between axisymmetric equilibria and nonaxisymmetric configurations? This is a challenging question. In the single SID case, it was suggested (Shu et al. 2000; Galli et al. 2001) that these solutions signal onsets of bifurcations from an axisymmetric SID to nonaxisymmetric SIDs (such as eccentric, oval, triangular distortions corresponding to $`m=1,2,3`$ etc.). Based on transmisions and overreflections of leading and/or trailing spiral density waves across the corotation in a time-dependent problem, Shu et al. (2000) also made a novel suggestion that the condition for stationary unaligned or logarithmic spiral configurations in a SID would determine whether spiral density waves may be swing-amplified (Goldreich & Lynden-Bell 1965; Fan & Lou 1997). In addition, Shu et al. (2000) examined the parameters of these stationarity conditions in reference to the disk stability criterion hypothesized by Ostriker & Peebles (1973) for the onset of bar-type instabilities (Miller et al. 1970; Hohl 1971; Kalnajs 1972) and estimated Ostriker-Peebles criterion for aligned secular and spiral dynamic barlike instabilities in a SID. Judging the similarity and difference between a single SID and a composite SID system, we examine stability properties for stationary aligned and unaligned configurations in the same spirit.

It should be noted that onsets of or transitions to barlike instabilities may also occur at non-zero pattern speeds such as those of Maclaurin spheroids to Jacobi ellipsoids (Chandrasekhar 1969) or those of Kalnajs disks to bar configurations etc (Binney & Tremaine 1987). The most relevant analogy here is perhaps the transitions from Maclaurin spheroids to stationary Dedekind ellipsoids with fixed configurations in space (e.g., Chandrasekhar 1969).

Let us first perform an analysis on aligned secular barlike instability in a composite SID system. Similar to the single SID case of Shu et al. (2000), the criterion postulated here is expressed in terms of the ratio of the kinetic energy of SID rotation $`𝒯`$ to the absolute value of the gravitational potential energy $`𝒲`$ (Ostriker & Peebles 1973; Binney & Tremaine 1987). In integral forms, we derive

$$𝒯_0^R\frac{1}{2}\mathrm{\Sigma }_0^s(r\mathrm{\Omega }_s)^22\pi r𝑑r+_0^R\frac{1}{2}\mathrm{\Sigma }_0^g(r\mathrm{\Omega }_g)^22\pi r𝑑r$$
(65)
and

$$𝒲_0^Rr(\mathrm{\Sigma }_0^s+\mathrm{\Sigma }_0^g)\frac{d\varphi _0}{dr}2\pi r𝑑r$$
(66)
by definitions, where $`R`$ is a radius that is allowed to go to infinity. To derive the virial theorem for the composite SID system from the background equilibrium conditions, we write equations (2) and (5) in the forms of

$$\mathrm{\Sigma }_0^sr\mathrm{\Omega }_s^2=\frac{d}{dr}(a_s^2\mathrm{\Sigma }_0^s)\mathrm{\Sigma }_0^s\frac{d\varphi _0}{dr}$$
(67)
and

$$\mathrm{\Sigma }_0^gr\mathrm{\Omega }_g^2=\frac{d}{dr}(a_g^2\mathrm{\Sigma }_0^g)\mathrm{\Sigma }_0^g\frac{d\varphi _0}{dr}.$$
(68)
Adding equations (67) and (68), we obtain

$$(\mathrm{\Sigma }_0^s\mathrm{\Omega }_s^2+\mathrm{\Sigma }_0^g\mathrm{\Omega }_g^2)r=\frac{d}{dr}(a_s^2\mathrm{\Sigma }_0^s+a_g^2\mathrm{\Sigma }_0^g)+(\mathrm{\Sigma }_0^s+\mathrm{\Sigma }_0^g)\frac{d\varphi _0}{dr}.$$
(69)
Multiplying equation (69) by $`2\pi r^2dr`$ and integrating from $`0`$ to $`R`$, we arrive at

$$2(𝒯+𝒰)+𝒲=2\pi R^2[a_s^2\mathrm{\Sigma }_0^s(R)+a_g^2\mathrm{\Sigma }_0^g(R)],$$
(70)
where

$$𝒰_0^R(a_s^2\mathrm{\Sigma }_0^s+a_g^2\mathrm{\Sigma }_0^g)2\pi r𝑑r$$
(71)
is the thermal energy contained in the composite SID system. Equation (70) stands for the virial theorem generalized for a composite SID system.

Using equilibrium equations $`(8)(11)`$, the two integrals (65) and (66) can be expressed as

$$𝒯=a_s^4(D_s^2+1)\frac{D_s^2+(D_s^2+11/\beta )\delta }{2G(1+\delta )}R$$
(72)
and

$$𝒲=\frac{a_s^4(D_s^2+1)^2}{G}R.$$
(73)
For a composite SID system of infinite radial extent, both integrals diverge as $`R\mathrm{}`$ but their ratio remains finite. So the ratio of the kinetic energy of rotation to the absolute value of the gravitational potential energy is

$$\begin{array}{c}\hfill \frac{𝒯}{|𝒲|}=\frac{D_s^2+(D_s^2+11/\beta )\delta }{2(1+D_s^2)(1+\delta )}=\frac{1}{2}\frac{1+\delta /\beta }{2(D_s^2+1)(1+\delta )}.\end{array}$$
(74)
Equations $`(72)(74)`$ can all be explicitly symmetrized with respect to the parameters of the two SIDs by using background condition (11). We here use equation (74) as an example to illustrate this symmetry between parameters of two SIDs, namely,

$$\begin{array}{c}\hfill \frac{𝒯}{|𝒲|}=\frac{a_s^2\mathrm{\Sigma }_sD_s^2+a_g^2\mathrm{\Sigma }_gD_g^2}{[a_s^2(1+D_s^2)+a_g^2(1+D_g^2)](\mathrm{\Sigma }_s+\mathrm{\Sigma }_g)}.\end{array}$$

Note that the value of $`𝒯/|𝒲|`$ falls between $`0`$ and $`1/2`$ as usual (Binney & Tremaine 1987) and increases with the increase of $`D_s^2`$. Therefore, for stationary aligned configurations in a composite system of two coupled SIDs, the two possible values of $`D_s^2`$ correspond to two values of $`𝒯/|𝒲|`$ ratio; the larger and smaller values of $`D_s^2`$ correspond to larger and smaller values of $`𝒯/|𝒲|`$ ratio.

To illustrate the physical significance of the above result, let us examine a specific case when $`|m|=2`$, $`\delta =0.25`$ and $`\beta =2`$. Equation (51) now yields two solutions $`D_s^2=1.0552`$ and $`D_s^2=0.0948`$ (see Fig. 1). Substitution of these two values of $`D_s^2`$ into equation (74) in order gives $`𝒯/|𝒲|=0.2810`$ and $`𝒯/|𝒲|=0.0890`$, respectively. Based on numerical simulation experiments (involving $`300`$body particles) for stability of a rotating disk, Ostriker & Peebles (1973) suggested empirically that the approximate condition $`𝒯/|𝒲|\stackrel{<}{}\mathrm{\hspace{0.33em}0.14}\pm 0.02`$ is necessary but not sufficient for stability against bar-type instabilities555Binney & Tremaine (1987) seems to suggest both necessity and sufficiency of this criterion.. That is, when $`𝒯/|𝒲|\stackrel{>}{}\mathrm{\hspace{0.33em}0.14}\pm 0.02`$, a disk system would rapidly evolve into bar-type configurations (Miller et al. 1970; Hohl 1971; Hunter 1977). By our numerical example above and the analogy to the single SID case (Shu et al. 2000), the transition criteria from an axisymmetric equilibrium to aligned secular barlike instabilities via two different modes (i.e., upper $`y_1`$ and lower $`y_2`$ solutions) correspond two different values of $`𝒯/|𝒲|`$ ratio (larger and smaller than $`0.14`$ respectively). These variations of $`𝒯/|𝒲|`$ ratio appear considerable but not totally surprising because for Kalnajs disks (Kalnajs 1972), the relevant criterion is $`𝒯/|𝒲|<0.1268`$ for the stability of bar modes (see also discussions of Binney & Tremaine 1987).

There are two possible interpretations for the $`y_2`$ solution with a fairly low $`𝒯/|𝒲|=0.0890`$. First, if the correspondence between the $`𝒯/|𝒲|`$ ratio and the onset of secular barlike instabilities holds as suggested by Shu et al. (2000) for a single SID case and if our extension of this correspondence to a composite SID system is valid, then our analysis suggests that the threshold of $`𝒯/|𝒲|`$ ratio can be lowered considerably in a composite SID system. Second, as emphasized by Ostriker & Peebles (1973) in their note added in proof, there can be other instabilities that occur for $`𝒯/|𝒲|\stackrel{<}{}\mathrm{\hspace{0.33em}0.14}\pm 0.02`$. If this is true, then the suspected correspondence between the $`𝒯/|𝒲|`$ ratio and the onset of secular barlike instabilities might be coincidental in the single SID case. We are inclined towards the first interpretation although numerical simulations involving a composite disk system are deemed necessary to resolve this important issue. At any rate, a composite SID system tends to be stabilized by introducing a sufficiently massive dark-matter halo for both solutions (see our analysis on composite partial SID system later).

For the same value of $`\beta =2`$ but with a smaller $`\delta `$ value, e.g., $`\delta =0.1`$, we then have $`D_s^2=0.8557`$ and $`D_s^2=0.2125`$, corresponding to $`𝒯/|𝒲|=0.2428`$ and $`𝒯/|𝒲|=0.1064`$, respectively. To go further for $`\beta =2`$ and $`\delta =0.05`$, we have $`D_s^2=0.7500`$ and $`D_s^2=0.2857`$, corresponding to $`𝒯/|𝒲|=0.2211`$ and $`𝒯/|𝒲|=0.1204`$, respectively. By this sequence of three numerical examples with fixed $`\beta =2`$ and $`m=2`$, the tendency is clear that the ratio $`𝒯/|𝒲|`$ for $`y_1`$ decreases but remains fairly high, while the ratio $`𝒯/|𝒲|`$ for $`y_2`$ gradually approaches the usual estimate of $`0.14`$ as $`\delta `$ becomes smaller.

We now complement these numerical examples by an analytical analysis. As defined by equation (74), the ratio of $`𝒯/|𝒲|`$ involves three parameters $`\delta `$, $`\beta `$ and $`D_s^2`$ explicitly, and increases with increasing values of $`D_s^2`$ and $`\beta `$. Once $`\delta `$ and $`\beta `$ are given, the marginal value of $`D_s^2`$ can be determined by equation (51) for the onset of aligned barlike instabilities with $`m=2`$. We investigate below variation trends of $`𝒯/|𝒲|`$ versus $`\delta `$ when $`\beta `$ is specified.

Let us first start with $`\beta =1.5`$. For $`m=2`$, the critical value $`\beta _c`$ as given by equation (56) becomes

$$\begin{array}{c}\hfill \beta _c=\frac{3}{2}\left[1+\frac{1}{4\delta +1}\right].\end{array}$$

The condition of $`\beta \beta _c`$ (see the lower branch $`y_2`$ solutions in Fig 1) requires that for a given $`\beta `$, $`\delta `$ must fall within a proper range so as to render the $`y_2`$ solution branch physically meaningful. It turns out that for $`\beta =1.5`$, $`\delta `$ can take arbitrary values from $`0`$ to $`\mathrm{}`$. We then solve equation (51) for aligned stationary perturbations with $`m=2`$ and $`\beta =1.5`$ to obtain the $`y_2=D_s^2`$ solution branch as an explicit function of $`\delta `$, namely

$$\begin{array}{c}\hfill y_2=\frac{3+4\delta (1+16\delta +16\delta ^2)^{1/2}}{4(1+\delta )}.\end{array}$$

Substitution of this $`y_2`$ into equation (74) with $`\beta =1.5`$ yields an expression for marginal ratio $`𝒯/|𝒲|`$ as a function of $`\delta `$, namely

$$\begin{array}{c}\hfill \frac{𝒯}{|𝒲|}=\frac{9+16\delta 3(1+16\delta +16\delta ^2)^{1/2}}{6[7+8\delta (1+16\delta +16\delta ^2)^{1/2}]},\end{array}$$

which has a minimum value of $`5/360.1389`$ at $`\delta =2/30.67`$ (see the upper curve in Fig. 2).

For a somewhat larger $`\beta `$ such as $`\beta =2`$, the value of $`\delta `$ is no longer arbitrary in order to obtain nonnegative $`y_2`$; it turns out that $`0<\delta 0.5`$. Repeating the same procedure for the $`\beta =1.5`$ case, we then derive another expression for critical ratio $`𝒯/|𝒲|`$ in terms of $`\delta `$ as

$$\begin{array}{c}\hfill \frac{𝒯}{|𝒲|}=\frac{411\delta +3(8\delta +9\delta ^2)^{1/2}}{6[45\delta +(8\delta +9\delta ^2)^{1/2}]}.\end{array}$$

Within the interval of $`0<\delta 0.5`$, this ratio attains a minimum value of $`1/120.0833`$ at $`\delta =0.5`$ (see the middle curve in Fig. 2).

By similar considerations for the case of $`\beta =2.5`$, the required range of $`\delta `$ becomes $`0<\delta <1/8`$, and the critical ratio $`𝒯/|𝒲|`$ attains the minimum value of 1/30 at $`\delta =1/8`$ (see the lower curve in Fig. 2).

With a further increase of $`\beta `$, the required range of $`\delta `$ diminishes. When $`\beta >3`$, there is no allowed value of $`\delta `$ that makes $`y_2`$ physically meaningful.

These three examples are shown in Figure 2 for the variations of the marginal ratio $`𝒯/|𝒲|`$ versus $`\delta `$ with different values of $`\beta `$ at a fixed value of $`m=2`$. Each curve corresponds to an allowed range of $`\delta `$ for $`y_2=D_s^2>0`$. Our analysis of aligned cases already indicate a more complicated stability properties of a composite SID system than those of a single SID system. In particular, possible modes of smallest $`D_s^2`$ can have different values of $`m`$ as shown in Figs. 3 and 4.

Starting from a specific example with $`\delta =0.2`$ and $`\beta =1.5`$, the $`y_2`$ solution branch of aligned equation (51) yields a curve of $`D_s^2`$ versus the azimuthal wavenumber $`m2`$ as shown in Fig. 3. Qualitatively, this curve differs from the result of one-SID case, where $`D^2=|m|/(|m|+2)`$ increases monotonically with increasing $`m2`$ (see equation (27) of Shu et al. 2000 and our equation (52)). Here, $`D_s^2`$ attains a maximum at $`m=2`$ and $`3`$ and a minimum at $`m=7`$. In this particular case, both $`m=2`$ and $`m=3`$ give the same value of $`D_s^2`$ according to aligned equation (51).

For the same $`\delta =0.2`$ but a smaller $`\beta =1.2`$, the curve of $`D_s^2`$ for $`y_2`$ solution branch of equation (51) versus $`m2`$ is shown in Fig. 4 with the smallest $`D_s^2`$ at $`m=2`$. Note the variation in solution structures in Figs. 3 and 4 for a slight change of $`\beta `$ value.

### 3.2 Unaligned or spiral disturbances

For stationary unaligned or spiral disturbances, we take the following set of exact density-potential relation,

$$\begin{array}{c}\hfill \mu ^s=\sigma ^sr^{3/2}\mathrm{exp}(i\alpha \mathrm{ln}r),\\ \hfill \mu ^g=\sigma ^gr^{3/2}\mathrm{exp}(i\alpha \mathrm{ln}r),\end{array}$$
(75)
$$V=\upsilon ^sr^{1/2}\mathrm{exp}(i\alpha \mathrm{ln}r)+\upsilon ^gr^{1/2}\mathrm{exp}(i\alpha \mathrm{ln}r),$$
(76)
where $`\sigma ^s`$, $`\sigma ^g`$, $`\alpha `$, $`\upsilon ^s`$, $`\upsilon ^g`$ are constant coefficients and $`\upsilon ^s`$ and $`\sigma ^s`$ and $`\upsilon ^g`$ and $`\sigma ^g`$ are related by

$$\upsilon ^s=2\pi G𝒩_m(\alpha )\sigma ^s,$$

$$\upsilon ^g=2\pi G𝒩_m(\alpha )\sigma ^g.$$

Here, $`𝒩_m(\alpha )K(\alpha ,m)`$ is the Kalnajs function (Kalnajs 1971). For spiral perturbations in a composite SID system, this seems to be a sensible extension of earlier work relevant to the subject (e.g., Lynden-Bell & Lemos 1993; Syer & Tremaine 1996; Shu et al. 2000; Lou 2002). Note that the radial scaling parameter $`\alpha `$ (related to the radial wavenumber) is naturally taken to be the same in both stellar and gaseous SIDs.

In our analysis and computations, we shall make use of two useful formula of $`𝒩_m(\alpha )`$. One is the recursion relation in $`m`$ of $`𝒩_m(\alpha )`$ for a fixed $`\alpha `$ (Kalnajs 1971),

$$𝒩_{m+1}(\alpha )𝒩_m(\alpha )=[(m+1/2)^2+\alpha ^2]^1$$
(77)
and the other is the asymptotic expression of $`𝒩_m(\alpha )`$ (Shu et al. 2000),

$$𝒩_m(\alpha )(m^2+\alpha ^2+1/4)^{1/2}$$
(78)
for $`m^2+\alpha ^21`$. When accuracy is not that crucial in some quantitative considerations, asymptotic expression (78) is also used to compute the $`|m|=1`$ spiral solutions later.

For $`|m|>0`$, we proceed to solve equations (39) and (40) which, with a little algebra, can be cast into the compact forms of

$$\mu ^s=\left(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2}\right)(H_1r\mu ^s+G_1r\mu ^g),$$
(79)
$$\mu ^g=\left(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2}\right)(H_2r\mu ^g+G_2r\mu ^s),$$
(80)
where coefficients $`H_1`$, $`H_2`$, $`G_1`$ and $`G_2`$ are defined by

$$\begin{array}{c}\hfill H_1\frac{1}{D_s^2(m^22)}\left[1\frac{(D_s^2+1)𝒩_m(\alpha )}{1+\delta }\right],\\ \hfill H_2\frac{1}{D_g^2(m^22)}\left[1\frac{(D_g^2+1)𝒩_m(\alpha )\delta }{1+\delta }\right],\end{array}$$
(81)
$$\begin{array}{c}\hfill G_1\frac{(D_s^2+1)}{D_s^2(m^22)}\frac{𝒩_m(\alpha )}{1+\delta },\\ \hfill G_2\frac{(D_g^2+1)}{D_g^2(m^22)}\frac{𝒩_m(\alpha )\delta }{1+\delta }.\end{array}$$
(82)
We substitute $`\mu ^s`$ and $`\mu ^g`$ in the forms of equation (75) and (76) into equations (79) and (80) to obtain

$$[1H_1(m^2+\alpha ^2+1/4)]\mu ^s=G_1(m^2+\alpha ^2+1/4)\mu ^g,$$
(83)
$$[1H_2(m^2+\alpha ^2+1/4)]\mu ^g=G_2(m^2+\alpha ^2+1/4)\mu ^s.$$
(84)
Without the gravitational coupling between the stellar and gaseous SIDs, the result from either disk would be the same as that of Shu et al. (2000) for a single SID.

By multiplying both sides of equations (83) and (84) and removing $`\mu ^s\mu ^g`$, we obtain

$`[1H_1(m^2+\alpha ^2+1/4)][1H_2(m^2+\alpha ^2+1/4)]`$
$`=G_1G_2(m^2+\alpha ^2+1/4)^2.`$ (85)
By equations (81) and (82), equation (85) becomes

$$\begin{array}{c}\hfill \left\{1\frac{1}{D_s^2(m^22)}\left[1\frac{(D_s^2+1)𝒩_m(\alpha )}{1+\delta }\right]\left(m^2+\alpha ^2+\frac{1}{4}\right)\right\}\\ \hfill \times \{1\frac{1}{D_g^2(m^22)}[1\frac{(D_g^2+1)𝒩_m(\alpha )\delta }{1+\delta }]\\ \hfill \times (m^2+\alpha ^2+\frac{1}{4})\}\\ \hfill =\frac{D_s^2+1}{D_s^2(m^22)}\frac{D_g^2+1}{D_g^2(m^22)}\frac{𝒩_m^2(\alpha )\delta }{(1+\delta )^2}\left(m^2+\alpha ^2+\frac{1}{4}\right)^2.\end{array}$$
(86)
Equation (86) can be rewritten in a physically more informative form of

$$\begin{array}{c}\hfill \{D_s^2(m^22)[1\frac{(D_s^2+1)𝒩_m(\alpha )}{1+\delta }]\\ \hfill \times (m^2+\alpha ^2+\frac{1}{4})\}\\ \hfill \times \{D_g^2(m^22)[1\frac{(D_g^2+1)𝒩_m(\alpha )\delta }{1+\delta }]\\ \hfill \times (m^2+\alpha ^2+\frac{1}{4})\}\\ \hfill =(D_s^2+1)(D_g^2+1)\frac{𝒩_m^2(\alpha )\delta }{(1+\delta )^2}\left(m^2+\alpha ^2+\frac{1}{4}\right)^2,\end{array}$$
(87)
which is an exact relation among $`D_s`$, $`D_g`$, $`m`$, $`\alpha `$ and $`\delta `$.

By substituting relation (11), namely $`D_g^2=\beta (D_s^2+1)1`$, into equation (87), we obtain a quadratic equation of $`yD_s^2`$,

$$C_2y^2+C_1y+C_0=0,$$
(88)
where three coefficients are defined by

$$\begin{array}{c}\hfill C_2=\beta (m^22)[m^22+(m^2+\alpha ^2+1/4)𝒩_m(\alpha )],\text{ }(88a)\end{array}$$

$$\begin{array}{c}\hfill C_1=[\frac{2+\delta }{1+\delta }(m^22)(m^2+\alpha ^2+\frac{1}{4})𝒩_m(\alpha )\frac{\delta }{1+\delta }\\ \hfill \times (m^2+\alpha ^2+\frac{1}{4})^2𝒩_m(\alpha )(m^22)(\alpha ^2+\frac{9}{4})]\beta (88b)\\ \hfill \left(2m^2+\alpha ^2\frac{7}{4}\right)\left[m^22+\frac{(m^2+\alpha ^2+1/4)𝒩_m(\alpha )}{1+\delta }\right],\text{ }\end{array}$$

$$\begin{array}{c}\hfill C_0=(m^2+\alpha ^2+\frac{1}{4})\{(m^22)[\frac{𝒩_m(\alpha )}{1+\delta }1]\frac{\delta }{1+\delta }\\ \hfill \times (m^2+\alpha ^2+\frac{1}{4})𝒩_m(\alpha )\}\beta +(m^2+\alpha ^2+\frac{1}{4})\\ \hfill \times (2m^2+\alpha ^2\frac{7}{4})[1\frac{𝒩_m(\alpha )}{1+\delta }].(88c)\end{array}$$

#### 3.2.1 The analytical case of $`\beta =1`$

Parallel to the analysis for aligned case, we first investigate the $`\beta =1`$ case to gain useful insight. When $`\beta =1`$ and thus $`D_g^2=D_s^2`$ by condition (11), equation (88) becomes

$$\begin{array}{c}\hfill [D_s^2(m^22)(m^2+\alpha ^2+1/4)]\\ \hfill \times \{D_s^2[(m^22)+𝒩_m(\alpha )(m^2+\alpha ^2+1/4)]\\ \hfill +[𝒩_m(\alpha )1](m^2+\alpha ^2+1/4)\}=0\end{array}$$
(89)
which has two branches of solution

$$D_s^2=y_1=\frac{m^2+\alpha ^2+1/4}{m^22}$$
(90)
and

$$D_s^2=y_2=\frac{[1𝒩_m(\alpha )](m^2+\alpha ^2+1/4)}{(m^22)+𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}.$$
(91)
Meanwhile from equation (83), we have

$`{\displaystyle \frac{\mu ^g}{\mu ^s}}={\displaystyle \frac{1H_1(m^2+\alpha ^2+1/4)}{G_1(m^2+\alpha ^2+1/4)}}`$
$`=1{\displaystyle \frac{[D_s^2(m^22)(m^2+\alpha ^2+1/4)](1+\delta )}{(D_s^2+1)𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}}.`$ (92)
For $`D_s^2=y_1`$ in equation (92), one gets $`\mu ^g/\mu ^s=1`$ and for $`D_s^2=y_2`$ in equation (92), one gets $`\mu ^g/\mu ^s=\delta `$. These results parallel those of the aligned case. Note that by solution (90), one has unphysical cases of $`D_s^2<0`$ for $`|m|=1`$. For the solution branch (91) with $`|m|=1`$, the situation is somewhat involved. We shall return to this case later at the end of subsection 3.2.3 around equations $`(100)(104)`$.

In general, $`\beta >1`$ and $`\delta `$ varies. For different values of $`\delta `$ and $`\beta `$ with $`|m|1`$, equation (88) or (89) always yields two solutions of $`D_s^2`$. For $`|m|=1`$, one solution is unphysical. For $`|m|2`$, we find that for the two solutions of $`D_s^2`$, the ratio of surface mass density perturbations $`\mu ^g/\mu ^s`$ is either positive or negative. This means that for stationary spiral disturbances in a composite SID system, $`\mu ^g`$ and $`\mu ^s`$ are either in phase or out of phase. We shall show details later.

Here we have just discussed the special analytic case of $`\beta =1`$ and point out that the $`|m|=1`$ case is unique. Our subsequent computation and analysis for general $`\beta `$ values are case-specific from the $`|m|=0`$ case to the $`|m|2`$ cases and then back to the $`|m|=1`$ case.

#### 3.2.2 Marginal stability of axisymmetric disturbances

While equation (88) or (89) is derived for cases of $`|m|>0`$, it can also be used to describe the marginal stability for the special case of $`|m|=0`$, that is, stationary axisymmetric disturbances. Equation (89) with $`|m|=0`$ yields two branches of solution constrained by the physical requirement666Not shown here is a second solution branch of $`D_s^2`$ which is negative and thus unphysical (see Figure $`D1`$ in Appendix D for an example corresponding to the same parameters of Figure 2 and also see equation (90) for $`\beta =1`$ with $`m=0`$). of $`D_s^20`$. For the positive portions of solution $`D_s^2`$ with not-so-extreme parameters, we find the basic $`D_s^2`$ versus $`\alpha `$ profile is qualitatively similar to Figure 2 of Shu et al. (2000) as expected. When $`D_s^2`$ falls in the range between the maximum of the collapse regime and the minimum of the ring fragmentation curve, the composite SID system can support stable oscillations as indicated in Figs $`59`$. In the limit of $`\alpha 0`$, stable “breathing modes” can be sustained (Lemos et al. 1991). When $`\delta \mathrm{\Sigma }_0^g/\mathrm{\Sigma }_0^s`$ is small, which means the gas mass is small compared to the stellar mass, the gas influence on the stellar disk is weak and the marginal stability will be slightly modified compared to that of a single SID situation. As $`\delta `$ increases, the role of gas disk becomes dynamically more important. For young disk galaxies during their early epochs of formation and evolution, the gas materials can be more abundant than or comparable to stellar mass, the axisymmetric marginal stability should be quite different from the single SID situation.

As $`\delta `$ and/or $`\beta `$ vary, the marginal stability curves change accordingly as shown in Figures $`510`$. For small $`\delta `$ and $`\beta 1`$, the shape of marginal stability curve is similar to that of a single stellar SID (Shu et al. 2000). The curve contains two separate branches in the $`D_s^2\alpha `$ plane. The first collapse branch starts at the vertical $`D_s^2>0`$ axis with $`\alpha =0`$ and goes down to the horizontal $`\alpha >0`$ axis with $`D_s^2=0`$. The second ring fragmentation branch has a vertical asymptote, where the value of $`D_s^2`$ approaches infinity at a finite $`\alpha `$. The curve descends with increasing $`\alpha `$ to a minimum and the ascends to infinity as $`\alpha \mathrm{}`$. There are no other positive values for $`D_s^2`$ that can be found between the two branches just described. When $`\delta `$ is fixed and $`\beta `$ takes larger values, which means the stellar velocity dispersion further exceeds the sound speed of the gas SID, both branches descend, but this trend of variation is more conspicuous for the ring fragmentation branch.

For $`\delta =0.2`$, the increase of $`\beta `$ significantly lowers the ring fragmentation curve. In comparison, the collapse branch changes slightly with increasing $`\beta `$ when $`\delta `$ is small. This is shown in Figs 5 and 6. As $`\delta `$ increases for fixed $`\beta `$ value, the collapse branch shrinks with the surrounded area reduced (see Figs. 6 and 7). The collapse regime continues to shrink with increasing $`\delta `$. When $`\delta `$ becomes sufficiently large, this branch will completely disappear from the $`D_s^2\alpha `$ diagram (see Figs. 8 and 9), leaving only the ring fragmentation branch (Fig. 9). Meanwhile for the ring fragmentation branch, the trend is clear that as $`\delta `$ increases, the minimum of the marginal stability curve decreases at larger values of $`\alpha `$. Eventually, the curve crosses the $`\alpha `$ axis and becomes negative when $`\delta `$ and $`\beta `$ take proper values (Fig 10). We conclude that as $`\delta `$ or $`\beta `$ increases, both marginal stability curves move downward (towards smaller $`D_s^2`$).

It is interesting to note that the location $`\alpha `$ of the vertical asymptote for the ring fragmentation branch remains fixed independent of $`\delta `$ and $`\beta `$. By equation (88), $`D_s^2`$ diverges when $`C_2`$ defined by equation $`(88a)`$ vanishes, or equivalently777This is also true when magnetic field is involved in a SID (Shu et al. 2000; Lou 2002; Lou & Fan 2002)., $`𝒩_0(\alpha )(\alpha ^2+1/4)=2`$. By recursion relation (77) and the asymptotic expression (78), we may choose approximately

$`𝒩_0(\alpha )=(9/4+\alpha ^2)/(1/4+\alpha ^2)𝒩_2(\alpha )`$
$`=(9/4+\alpha ^2)/[(1/4+\alpha ^2)(17/4+\alpha ^2)^{1/2}]`$ (93)
to estimate an $`\alpha `$ value of $`1.7928`$ for $`C_2=0`$ in equation (88) (Shu et al. 2000; Lou 2002; Lou & Fan 2002).

We now provide physical interpretations for the two curves of marginal stability. The parameter $`\alpha `$ is a measure of radial wavenumber. Small $`\alpha `$ corresponds to larger radial perturbation scale. Therefore the composite SID system is vulnerable to Jeans’ instability when $`\alpha `$ is sufficiently small. Because of the angular momentum conservation, a SID rotation tends to work against gravitational collapse. Therefore, it requires a smaller $`\alpha `$ to induce Jeans’ collapse in the presence of disk rotation. When disk rotation is sufficiently fast ($`D_s^2`$ larger than a critical value at $`\alpha =0`$), Jeans’ collapse can be completely prevented.

On the other hand, the composite SID system should be also vulnerable to axisymmetric Toomre-type instability.888The original criterion was derived for local axisymmetric perturbations in a rotating disk (Safronov 1960; Toomre 1964). In various galactic contexts, it is of considerable interest to search for an effective $`Q`$ parameter for a composite disk system (Jog & Solomon 1984a,b; Bertin & Romeo 1988; Kennicutt 1989; Romeo 1992; Wang & Silk 1994; Elmegreen 1995; Jog 1996; Lou & Fan 1998b). Shu et al. (2000) noticed that the minimum of the ring fragmentation branch is effectively related to the Toomre $`Q`$ parameter (Toomre 1964) in a very close manner. In a magnetized SID with a coplanar magnetic field, Lou (2002) and Lou & Fan (2002) noticed that that the minimum of the MHD ring fragmentation branch is effectively related to the generalized MHD $`Q_M`$ parameter (Lou & Fan 1998a) in a very close manner. Physically, we therefore expect that the minimum of the ring fragmentation branch should be intimately related to an effective $`Q`$ parameter in a composite disk system (Elmegreen 1995; Jog 1996; Lou & Fan 1998b). When $`D_s^2`$ exceeds the minimum of the ring fragmentation branch, the composite system of two coupled SIDs becomes unstable to ringlike fragmentation for the range of radial wavenumbers $`\alpha `$ between the two points of intersection of the horizontal line and the the marginal stability curve.

#### 3.2.3 Stationary logarithmic spiral configurations

Following the same procedure of obtaining solutions for axisymmetric disturbances, we derive solution curves for cases of $`|m|=1,2,3\mathrm{}`$. As expected, for each set of $`\{m`$, $`\delta `$, $`\beta \}`$, there are two solution curves for $`D_s^2`$ versus $`\alpha `$ in a composite SID system according to equation (88) or (89). In some cases, part of the lower solution may become negative which is unphysical. Figures $`1113`$ show a set of numerical examples with $`|m|=2`$.

In the case of $`|m|=1`$, only one branch of $`D_s^2`$ is physically meaningful. We will discuss more on this problem later. Let us first examine cases of $`|m|2`$.

For given values of $`\delta `$ and $`\beta `$, there are two solution curves in general for $`D_s^2`$ versus $`\alpha `$ when $`|m|2`$. The lower branch may cross the horizontal $`\alpha `$ axis when $`\delta `$ and/or $`\beta `$ are sufficiently large. In contrast, the upper branch will never become negative. We note that variations of $`\beta `$ and $`\delta `$ play distinctly different roles. The increase of $`\delta `$ mainly moves the upper branch upward and the lower branch downward, whereas the increase of $`\beta `$ tends to move both branches downward. For very small $`\delta `$ values, different $`\beta `$ values may lead to quite different solution configurations (see Figs. $`1113`$).

In order to fully understand solution behaviors (Figs. $`1113`$), we perform analytic analysis on properties of stationary spiral configurations as in the aligned cases of $`|m|2`$. Firstly, through extensive numerical experiments, one can show that $`yD_s^2`$ increases with wavenumber $`\alpha `$ for both upper and lower branches. Secondly, once $`|m|`$ and $`\delta `$ are known with fixed $`\alpha `$, we would like to know how $`yD_s^2`$ and $`\mu ^g/\mu ^s`$ vary with $`\beta `$. It turns out that for both solution branches, $`yD_s^2`$ monotonically decrease with increasing $`\beta `$ and $`\mu ^g/\mu ^s`$ decreases with increasing $`D_s^2`$. It is then clear that $`\mu ^g/\mu ^s`$ goes with $`\beta `$ in the same trend. In more informative forms, we derive two useful inequalities for $`y_1`$ and $`y_2`$ below. The first one for $`y_1`$ is thus

$$\begin{array}{c}\hfill (m^2+\alpha ^2+1/4)/(m^22)\\ \hfill \times \left\{1\frac{[m^22+(m^2+\alpha ^2+1/4)]𝒩_m(\alpha )}{[m^22+(m^2+\alpha ^2+1/4)𝒩_m(\alpha )](1+\delta )}\right\}\\ \hfill <y_1<(m^2+\alpha ^2+1/4)/(m^22)\end{array}$$
(94)
with the lower bound on the left-hand side determined by taking $`\beta \mathrm{}`$ and the upper bound on the right-hand side obtained by taking $`\beta =1`$; and the second one for $`y_2`$ is

$$1<y_2<\frac{[1𝒩_m(\alpha )](m^2+\alpha ^2+1/4)}{m^22+𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}$$
(95)
with the lower bound on the left-hand side determined by taking $`\beta \mathrm{}`$ and the upper bound on the right-hand side obtained by taking $`\beta =1`$. One can readily show that the lower bound of $`y_1`$ is positive definite as $`𝒩_m(\alpha )<1`$ for $`|m|2`$. Therefore, the $`D_s^2y_1`$ branch is always physically valid given values of $`\delta `$ and $`\beta `$ with $`|m|2`$. Moreover, as in the aligned cases, the left-hand side of inequality (94) for $`y_1`$ remains always larger than the right-hand side of inequality (95) for $`y_2`$ such that the two solution branches $`y_1`$ and $`y_2`$ will never intersect with each other because $`y_1>y_2`$ always. Meanwhile, there exists also a critical value $`\beta _c`$ for $`y_2D_s^2`$ branch at a given wavenumber $`\alpha `$ such that $`y_2D_s^2`$ becomes zero. In analytic form, this critical value $`\beta _c`$ is given by

$$\begin{array}{c}\hfill \beta _c=\frac{m^22+(m^2+\alpha ^2+1/4)}{m^22+𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}\{1+\\ \hfill \frac{𝒩_m(\alpha )[1𝒩_m(\alpha )](m^2+\alpha ^2+1/4)}{𝒩_m(\alpha )(m^2+\alpha ^2+1/4)\delta +(m^22)[1+\delta 𝒩_m(\alpha )]}\}.\end{array}$$
(96)
Since $`𝒩_m(\alpha )<1`$ for $`|m|2`$, the critical value $`\beta _c`$ remains always greater than $`1`$ and decreases with increasing $`\delta `$. When $`\delta \mathrm{}`$, the critical value $`\beta _c`$ approaches a limit of

$$\beta _{cLim}\frac{m^22+(m^2+\alpha ^2+1/4)}{m^22+𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}.$$
(97)
In parallel to the study of aligned cases, we further derive the range of $`\mu ^g/\mu ^s`$ for both solution branches when $`|m|`$ and $`\delta `$ are known at a given wavenumber $`\alpha `$. For the $`y_1D_s^2`$ branch, we thus have

$$1<\frac{\mu ^g}{\mu ^s}<\frac{(m^2+\alpha ^2+1/4)𝒩_m(\alpha )\delta }{(m^2+\alpha ^2+1/4)𝒩_m(\alpha )\delta +(m^22)(1+\delta )}$$
(98)
with the lower bound on the left-hand side determined by $`\beta =1`$ and the upper bound on the right-hand side determined by $`\beta \mathrm{}`$; for the $`y_2D_s^2`$ solution branch, we have

$$\delta <\frac{\mu ^g}{\mu ^s}<\frac{1+\delta }{𝒩_m(\alpha )}1$$
(99)
with the lower bound on the left-hand side determined by $`\beta =1`$ and the upper bound on the right-hand side determined by the critical value $`\beta _c`$ such that $`y_2D_s^2=0`$. Since $`\mu ^g/\mu ^s`$ increases with increasing $`\beta `$, we always have negative $`\mu ^g/\mu ^s`$ for the $`y_1D_s^2`$ branch and positive $`\mu ^g/\mu ^s`$ for the $`y_2D_s^2`$ branch.

It is remarkable that the spiral cases are strikingly similar to the aligned cases. If one replaces $`𝒩_m(\alpha )`$ by $`(m^2+\alpha ^2+1/4)^{1/2}`$ and $`(m^2+\alpha ^2+1/4)^{1/2}`$ by $`|m|`$, all the above results for the spiral cases will degenerate to the aligned cases. By these considerations of aligned and spiral cases with $`|m|2`$, this apparent similarity should not be surprising because the common physical nature of stationary aligned and spiral disturbances, that is, purely azimuthal and general density waves in a composite SID system balanced by disk rotation.

As promised earlier, we now get back to the $`|m|=1`$ case which is special for a stationary spiral configuration as well as for a stationary aligned configuration. For aligned $`|m|=1`$ case, equation (51) is satisfied for arbitrary $`D_s^2>0`$ (see Shu et al. 2000 for a single SID). For spiral $`|m|=1`$ case, when $`\delta `$ and $`\alpha `$ are given, the $`y_1`$ and $`y_2`$ solution branches will increase and decrease with increasing $`\beta `$, respectively. When $`\beta \mathrm{}`$, the two solutions are

$$y=1$$

and

$$y=\left(\alpha ^2+\frac{5}{4}\right)+\frac{𝒩_1(\alpha )(\alpha ^2+5/4)(\alpha ^2+1/4)}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}.$$
(100)
Note in the above equation, we have $`𝒩_1(\alpha )(\alpha ^2+5/4)1`$ always positive according to equations (77) and (78).

The determination of the specific bounds of these two solutions will depend on a relation between $`\delta `$ and $`\alpha `$. When

$$\begin{array}{c}\hfill y=\left(\alpha ^2+\frac{5}{4}\right)+\frac{𝒩_1(\alpha )(\alpha ^2+5/4)(\alpha ^2+1/4)}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}\\ \hfill =1\left(\alpha ^2+\frac{1}{4}\right)\frac{\delta [𝒩_1(\alpha )(\alpha ^2+5/4)1]1}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}>1\end{array}$$

which means $`\delta <[𝒩_1(\alpha )(\alpha ^2+5/4)1]^1`$, we have

$$(\alpha ^2+5/4)<y_1<1$$
(101)
and

$$\begin{array}{c}\hfill \left(\alpha ^2+\frac{5}{4}\right)+\frac{𝒩_1(\alpha )(\alpha ^2+5/4)(\alpha ^2+1/4)}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}\\ \hfill <y_2<\frac{[1𝒩_1(\alpha )](\alpha ^2+5/4)}{𝒩_1(\alpha )(\alpha ^2+5/4)1}.\end{array}$$
(102)
On the other hand, when

$$\begin{array}{c}\hfill y=\left(\alpha ^2+\frac{5}{4}\right)+\frac{𝒩_1(\alpha )(\alpha ^2+5/4)(\alpha ^2+1/4)}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}\\ \hfill =1\left(\alpha ^2+\frac{1}{4}\right)\frac{\delta [𝒩_1(\alpha )(\alpha ^2+5/4)1]1}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}<1\end{array}$$

which means $`\delta >[𝒩_1(\alpha )(\alpha ^2+5/4)1]^1`$, we have

$$\begin{array}{c}\hfill \alpha ^2\frac{5}{4}<y_1<\alpha ^2\frac{5}{4}+\frac{𝒩_1(\alpha )(\alpha ^2+5/4)(\alpha ^2+1/4)}{[𝒩_1(\alpha )(\alpha ^2+5/4)1](1+\delta )}\end{array}$$
(103)
and

$$1<y_2<\frac{[1𝒩_1(\alpha )](\alpha ^2+5/4)}{𝒩_1(\alpha )(\alpha ^2+5/4)1}.$$
(104)
Note that for bounds of these inequalities, one side is determined by $`\beta =1`$ and the other side by $`\beta \mathrm{}`$ for both solution branches. It becomes clear that spiral $`|m|=1`$ case differs from spiral $`|m|2`$ cases. In particular, the ranges of $`y_1`$ and $`y_2`$ solutions are controlled by a relation between $`\delta `$ and $`\alpha `$, that is, whether $`\delta >[𝒩_1(\alpha )(\alpha ^2+5/4)1]^1`$ or not, such that the two solution branches do not intersect with each other. By inequalities (101) and (103), the $`y_1`$ solution branch is always negative and thus unphysical. In contrast, by inequalities (102) and (104), it is still possible for portions of the $`y_2`$ solution branch to be positive when $`\delta `$, $`\alpha `$ and $`\beta `$ take on proper values. Note that this $`y_2`$ solution branch is the counterpart of a single SID case.

As indicated in Figs. $`1113`$, for given values of $`|m|`$, $`\beta `$ and $`\delta `$, it possible for a composite SID system to sustain upper (smaller $`\alpha `$) and lower (larger $`\alpha `$) solution branches with different values of $`\alpha `$ simultaneously. Likewise, one may switch the role of $`|m|`$ and $`\alpha `$, that is, given $`\alpha `$, a composite SID system can support various spiral solutions with different values of $`|m|`$ (see Fig. 3 of Shu et al. 2000). This is to be understood that both $`|m|`$ and $`\alpha `$ are characteristic properties of density waves rather than background parameters of a composite SID system. In the linear regime, this allowed superposition of different solutions can give rise to a variety of galactic pattern structures that cannot fit by one simple spiral pattern.

The interpretation for the stability properties of these stationary logarithmic spiral configurations is not so straightforward even for a single SID case because of the power-law background divergence at $`r0`$ and thus the difficulty of imposing proper boundary conditions there (Zang 1976; Toomre 1977; Lynden-Bell & Lemos 1999; Goodman & Evans 1999; Shu et al. 2000; Galli et al. 2001). Shu et al. (2000) speculated that the stationary spiral solution condition represents onset of spiral instabilities on the ground that transmission and overreflection of leading/trailing spiral density waves across the corotation can be swing-amplified. This was done mainly through numerical experiments and appeared to be consistent with an earlier proposed criterion (Goodman & Evans 1999). Furthermore, the stability of stationary spiral configurations may be examined in light of the Ostriker-Peebles criterion (Shu et al. 2000). It is then plausible to extend this interpretation to a composite SID system where apparently more leading/trailing spiral density waves across the corotation would be involved in the processes of transmission and overreflection.

Therefore, in Figs. $`1113`$, as $`D_s^2`$ increases from $`0`$ and first intersects with the lower solution branch, the lower branch spiral solution can be sustained and be vulnerable to dynamical barred-spiral instabilities. As $`D_s^2`$ increases further to intersect the upper solution branch, the upper branch spiral solution (together with the lower branch solution at a larger $`\alpha `$) can be sustained and be vulnerable to spiral instabilities. In the case of Fig. 11, $`D_s^2`$ may be small enough to avoid dynamical barred-spiral instabilities. But in cases of Figs. 12 and 13, dynamical barred-spiral instabilities just cannot be avoided for $`D_s^20`$ when $`\alpha `$ is sufficiently large.

For the specific case of Fig. 11 with $`m=2`$, $`\delta =0.2`$ and $`\beta =1.5`$, the two minima of the $`y_1`$ and $`y_2`$ solution branches, i.e. the values of $`D_s^2`$ when $`\alpha 0`$, are 1.3406 and 0.3697, respectively. The corresponding ratio of $`𝒯/|𝒲|`$ defined by equation (74) are 0.2982 and 0.1552 in order. In analogy to the single SID case of Shu et al. (2000), the two minima of $`D_s^2`$ signal onsets of dynamical barred-spiral instability in the composite SID system. By this interpretation, when the rotation of the stellar SID $`D_s^2`$ is smaller than 0.3697, the composite SID system can be stable against $`m=2`$ spiral instabilities. Otherwise, dynamical barred-spiral instabilities will develop in the composite SID system, similar to processes shown in sections 5.7 and 7 of Shu et al. (2000). Notice the introduction of the gaseous disk has decreased somewhat the marginal value of $`D_s^2`$ and thus the threshold ratio of $`𝒯/|𝒲|`$ in comparison to a single SID case.

Similar to the aligned cases, spiral stability properties of a composite SID system are more complicated than those of a single SID. For stationary logarithmic spiral configurations, $`D_s^2`$ usually attains a minimum value at $`\alpha =0`$ for both upper and lower solution branches of spiral equation (88), except for cases when $`y_2=D_s^2`$ solutions become negative for small $`\alpha `$ values (see Figs. 5-9).

As in the aligned case, we start from the same parameter set of $`\delta =0.2`$ and $`\beta =1.5`$ yet with a fixed $`\alpha =0`$. The $`y_2`$ solution branch of spiral equation (88) then yields a curve of $`D_s^2`$ versus $`m2`$ as shown in Fig. 14 where the smallest value of $`D_s^2=y_2`$ occurs at azimuthal periodicity $`m=9`$.

For the same $`\delta =0.2`$ and $`\alpha =0`$ yet a smaller $`\beta =1.2`$, the curve of $`D_s^2`$ for $`y_2`$ solution branch of spiral equation (88) versus $`m2`$ is shown in Fig. 15, where $`m=2`$ gives the smallest $`D_s^2`$.

## 4 Composite partial SID system

In this section, we consider stationary perturbations in a composite system of two gravitationally coupled partial SIDs (Syer & Tremaine 1996; Shu et al. 2000; Lou 2002; Lou & Fan 2002) to model the gravitational influence of a background axisymmetric dark-matter halo on the composite SID system. The response of this massive dark-matter halo to perturbations in the composite SID system is ignored. Formally, we introduce an additive gravity term $`\mathrm{\Phi }/r`$ or $`\mathrm{\Phi }/\phi `$ to the basic equations (2), (3) and (5), (6), where $`\mathrm{\Phi }`$ represents the gravitational potential contribution from the background dark-matter halo.

For momentum equations in the stellar disk, we have

$$\frac{u^s}{t}+u^s\frac{u^s}{r}+\frac{j^s}{r^2}\frac{u^s}{\phi }\frac{j^{s2}}{r^3}=\frac{1}{\mathrm{\Sigma }^s}\frac{}{r}(a_s^2\mathrm{\Sigma }^s)\frac{\varphi }{r}\frac{\mathrm{\Phi }}{r}$$
(105)
and

$$\frac{j^s}{t}+u^s\frac{j^s}{r}+\frac{j^s}{r^2}\frac{j^s}{\phi }=\frac{1}{\mathrm{\Sigma }^s}\frac{}{\phi }(a_s^2\mathrm{\Sigma }^s)\frac{\varphi }{\phi }\frac{\mathrm{\Phi }}{\phi }.$$
(106)
Likewise, for momentum equations in the gas disk, we have

$$\frac{u^g}{t}+u^g\frac{u^g}{r}+\frac{j^g}{r^2}\frac{u^g}{\phi }\frac{j^{g2}}{r^3}=\frac{1}{\mathrm{\Sigma }^g}\frac{}{r}(a_g^2\mathrm{\Sigma }^g)\frac{\varphi }{r}\frac{\mathrm{\Phi }}{r}$$
(107)
and

$$\frac{j^g}{t}+u^g\frac{j^g}{r}+\frac{j^g}{r^2}\frac{j^g}{\phi }=\frac{1}{\mathrm{\Sigma }^g}\frac{}{\phi }(a_g^2\mathrm{\Sigma }^g)\frac{\varphi }{\phi }\frac{\mathrm{\Phi }}{\phi }.$$
(108)
For a composite partial SID system, we now introduce a dimensionless $``$ parameter defined by $`\varphi /(\varphi +\mathrm{\Phi })`$. Now the background equilibrium should be modified accordingly. As before, we write $`\mathrm{\Omega }_s=a_sD_s/r`$, $`\mathrm{\Omega }_g=a_gD_g/r`$, and thus $`\kappa _s=\sqrt{2}\mathrm{\Omega }_s`$, $`\kappa _g=\sqrt{2}\mathrm{\Omega }_g`$. The equilibrium surface mass densities now become

$$\begin{array}{c}\hfill \mathrm{\Sigma }_0^s=\frac{a_s^2}{2\pi Gr}\frac{(1+D_s^2)}{1+\delta },\\ \hfill \mathrm{\Sigma }{}_{0}{}^{g}=\frac{a_g^2}{2\pi Gr}\frac{(1+D_g^2)\delta }{1+\delta }.\end{array}$$
(109)
Note that condition (11), namely, $`a_s^2(D_s^2+1)=a_g^2(D_g^2+1)`$ still holds. We refer to partial SIDs for $`0<<1`$, in contrast to full SIDs with $`=1`$.

The linearized coplanar perturbation equations should take the same forms of those in full SIDs written down in section 2.3. We thus follow the procedure of the full SIDs case up to equations (37) and (38). For the modified equilibrium with $`\omega =0`$, we have for the stellar partial SID

$$\begin{array}{c}\hfill m[\mu ^s+\frac{1}{D_s^2(m^22)}(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2})\\ \hfill \times (r\mu ^s+\frac{1+D_s^2}{2\pi G}\frac{V}{1+\delta })]=0,\end{array}$$
(110)
and for the gaseous SID

$$\begin{array}{c}\hfill m[\mu ^g+\frac{1}{D_g^2(m^22)}(\frac{m^2}{r}2\frac{d}{dr}r\frac{d^2}{dr^2})\\ \hfill \times (r\mu ^g+\frac{1+D_g^2}{2\pi G}\frac{\delta V}{1+\delta })]=0\end{array}$$
(111)
(see equations (39) and (40)). Equations (110) and (111) are to be solved simultaneously with Poisson’s integral (32).

We now consider consequences of a composite partial SID system in comparison with a composite system of two full SIDs.

For aligned cases: the $`|m|=0`$ disturbance is still a rescaling of one axisymmetric equilibrium to a neighboring axisymmetric equilibrium. For aligned cases of $`|m|=1`$, the generalized counterpart of equation (51) yields only one degenerate solution $`y=1`$ which is unphysical. In other words, eccentric displacements of $`|m|=1`$ are not allowed in a composite partial SID system. With $`\delta `$ and $`|m|`$ given, $`|m|2`$ cases yield two solutions satisfying the following inequalities

$$\begin{array}{c}\hfill \frac{m^2}{m^22}\frac{2|m|(m^21)}{(m^22)(m^2+|m|2)(1+\delta )}\\ \hfill <y_1<m^2/(m^22)\end{array}$$
(112)
and

$$1<y_2<|m|(|m|)/(m^2+|m|2),$$
(113)
respectively (see inequalities (54) and (55) for the cases of full SIDs). For a composite system of two coupled full SIDs with $`=1`$, these results naturally reduce to the familiar ones, namely, inequalities (54) and (55), derived earlier.

Given parameters $`\delta `$ and $`\beta `$ with $`|m|2`$, partial SIDs lead to two solution branches of $`D_s^2`$ that are larger than those of full SIDs. Therefore, the critical value $`\beta _c`$ for $`y_2D_s^2=0`$ is larger in partial SIDs, namely

$$\begin{array}{c}\hfill \beta _c=2(m^21)/(m^2+|m|2)\\ \hfill \times \left[1+\frac{|m|(|m|1)}{|m|(m^2+|m|2)\delta +(m^22)(|m|)}\right].\end{array}$$
(114)
For full SIDs with $`=1`$, expression (114) reduces to expression (56) as expected. Meanwhile, the ratio of surface mass density perturbations $`\mu ^g/\mu ^s`$ becomes

$$\frac{\mu ^g}{\mu ^s}=1\frac{[D_s^2(m^22)m^2](1+\delta )}{|m|(D_s^2+1)}.$$
(115)
Specifically for the two solution branches of $`D_s^2`$, we derive two inequalities, namely

$$1<\frac{\mu ^g}{\mu ^s}<\frac{|m|\delta }{(m^2+|m|2)\delta +m^22}$$
(116)
for the $`D_s^2=y_1`$ solution branch, and

$$\delta <\frac{\mu ^g}{\mu ^s}<|m|(1+\delta )/1$$
(117)
for the $`D_s^2=y_2`$ branch. Notice that expressions $`(115)`$, $`(116)`$ and $`(117)`$ have their obvious counterparts of the full SIDs cases of expressions $`(53)`$, $`(58)`$ and $`(59)`$ derived earlier.

The same procedure can be applied to the unaligned situation, which is expected to yield the same trends as the aligned situation.

We now derive the virial theorem in a composite partial SID system. By equations (105) and (107) for the equilibrium state, we obtain

$$\mathrm{\Sigma }_0^sr\mathrm{\Omega }_s^2=\frac{d}{dr}(a_s^2\mathrm{\Sigma }_0^s)\mathrm{\Sigma }_0^s\frac{d(\varphi _0+\mathrm{\Phi })}{dr}$$
(118)
and

$$\mathrm{\Sigma }_0^gr\mathrm{\Omega }_g^2=\frac{d}{dr}(a_g^2\mathrm{\Sigma }_0^g)\mathrm{\Sigma }_0^g\frac{d(\varphi _0+\mathrm{\Phi })}{dr}.$$
(119)
Adding equations $`(118)`$ and $`(119)`$, multiplying the resulting equation by $`2\pi r^2dr`$ and integrating from $`0`$ to $`R`$, we obtain

$$2(𝒯+𝒰)+𝒲=2\pi R^2[a_s^2\mathrm{\Sigma }_0^s(R)+a_g^2\mathrm{\Sigma }_0^g(R)],$$
(120)
where

$$𝒯_0^R\frac{1}{2}\mathrm{\Sigma }_0^s(r\mathrm{\Omega }_s)^22\pi r𝑑r+_0^R\frac{1}{2}\mathrm{\Sigma }_0^g(r\mathrm{\Omega }_g)^22\pi r𝑑r,$$
(121)
$$𝒰_0^R(a_s^2\mathrm{\Sigma }_0^s+a_g^2\mathrm{\Sigma }_0^g)2\pi r𝑑r,$$
(122)
$$𝒲_0^Rr(\mathrm{\Sigma }_0^s+\mathrm{\Sigma }_0^g)\frac{d(\varphi _0+\mathrm{\Phi })}{dr}2\pi r𝑑r$$
(123)
are the kinetic energy of rotation, the thermal energy and the gravitational potential energy, respectively. Note here the dark-matter halo contributes to the gravitational potential energy and $`\varphi _0/(\varphi _0+\mathrm{\Phi })`$.

## 5 Summary and discussions

There are four main tracks of research that lead to our investigation here on stationary perturbation configurations in a composite system of two gravitationally coupled full or partial SIDs. The first track is the classical study on spheroidal and ellipsoidal figures of equilibrium (Chandrasekhar 1969; Tassoul 1978) named after famous mathematicians including Maclaurin, Jacobi, Dedekind, Riemann and Poincaré. These spheroids and ellipsoids have their analogous compressible counterparts (Ostriker 1978), disk counterparts by collapsing along the rotation axis (Weinberg & Tremaine 1983; Weinberg 1983), as well as counterparts of composite systems (Vandervoort 1991a,b). The second track is along the past investigation of possible SID configurations and their stability properties (Zang 1976; Toomre 1977; Lemos, Kalnajs & Lynden-Bell 1991; Lynden-Bell & Lemos 1993; Lee & Goodman 1999; Evans & Read 1999; Goodman & Evans 1999; Syer & Tremaine 1996; Shu et al. 2000; Galli et al. 2001). The third track follows the studies of density waves in a composite disk system (Lin & Shu 1966; Jog & Solomon 1984a,b; Bertin & Romeo 1988; Elmegreen 1995; Jog 1996; Lou & Fan 1998b, 2000a,b). The fourth track pursues a basic understanding of MHD density waves in magnetized disk systems (Fan & Lou 1996; Lou & Fan 1998a, 2002, 2003; Shu et al. 2000; Lou, Yuan & Fan 2001; Lou 2002). Naturally, we shall further study possible stationary configurations in a composite system of gravitationally coupled stellar SID with gaseous MSID. These studies are of interests in galactic dynamics (Bertin & Lin 1996), star formation (Shu et al. 1999), the light cusps seen in the nuclei of galaxies (Crane et al. 1993), circumnuclear starburst “rings” and disk accretion processes around active galactic nuclei (Lou et al. 2001).

In reference to the work of Lou & Fan (1998b) on density waves in a composite disk system, of Shu et al. (2000) on stationary perturbation configurations of a single isopedically magnetized SID, and of Lou (2002) on stationary logarithmic fast and slow configurations in a single MSID, we studied here possible stationary configurations in a composite system of two gravitationally coupled full or partial SIDs using the two-fluid formalism, solved the Poisson equation exactly, and derived both aligned and unaligned or spiral stationary perturbation solutions. We have reached several conclusions and suggestions summarized below.

1. Aligned Stationary Configurations

For perturbation configurations in a composite SID system, there exist in general two branches of solutions as expected (Lou & Fan 1998b). This is one major difference between a composite SID system and a single SID case. We shall start with the aligned cases first.

The aligned axisymmetric $`|m|=0`$ case is again a rescaling from one axisymmetric equilibrium to a neighboring axisymmetric equilibrium (Shu et al. 2000; Lou 2002), except that now the rescaling occurs simultaneously in both SIDs.

We find as in the single full SID case (Shu et al. 2000) that, for eccentric $`|m|=1`$ displacements, equation (51) can be satisfied for unconstrained positive values of $`D_s^2`$. Given known idealizations of the full SID model, we consider that such eccentric displacements are possible. However, the situation is quite different for a composite system of two partial SIDs for which eccentric displacements are not allowed consistent with the similar result of a single partial SID (Lou 2002).

For $`|m|2`$, we derive two branches of $`D_s^2`$ solution for possible values of rotation speed parameter $`D_s`$ such that aligned stationary configurations are sustained in a composite SID system. Of the two branches, one is always larger than the other. For the larger branch denoted by $`y_1`$, perturbation enhancements of surface mass densities in stellar and gaseous disks are out of phase, while for the smaller branch denoted by $`y_2`$, the two enhancements are in phase. By the physical requirement of $`D_s^2>0`$, the larger $`y_1`$ branch is always valid (see inequality 54), while the smaller $`y_2`$ branch may become negative and thus unphysical (see inequality 55). Given $`|m|2`$ and $`\delta `$, when $`\beta `$ exceeds a critical value $`\beta _c`$ defined by equation (56), the smaller $`y_2`$ branch turns negative (see Fig. 1). In other words, only for disk galaxies with $`\beta <\beta _c`$, the smaller $`y_2`$ branch can be physically meaningful. Moreover, the smaller $`y_2`$ branch corresponds to the solution of single SID case as it would approach the limit of the single SID result when $`\beta 1`$, that is, when the system can be effectively viewed as a single SID. In contrast, the larger $`y_1`$ branch obtained in a composite SID system has no counterpart in the single SID case. Physically, stationary aligned disturbances actually involve purely azimuthal propagation of density waves counter balanced by SID rotation (Lou 2002; Lou & Fan 2002). The same density-wave interpretation holds for stationary unaligned or spiral perturbations. As in the single SID case, stationary aligned configurations here mark onsets of secular instabilities from an axisymmetric equilibrium to nonaxisymmetric configurations. For different values of $`|m|`$, distinct aligned configurations may coexist.
2. Unaligned or Spiral Stationary Configurations

Solution properties of the unaligned or spiral logarithmic configurations bear strong resemblance to those of the aligned cases. The spiral case involves both radial and azimuthal density wave propagations, whereas the aligned case involves only an azimuthal propagation of density waves. To a certain extent, we may regard aligned perturbation cases as special examples of spiral perturbation cases with $`|m|2`$.

For the marginal case of $`|m|=0`$, we have computed the marginal stability curves of $`D_s^2`$ versus radial wavenumber $`\alpha `$, invoking the same physical interpretations of Shu et al. (2000) for regimes of collapse, stable oscillations, and ring fragmentation. We shall further demonstrate the correctness of this interpretation by performing a WKBJ analysis in a separate paper. For different parameters $`\delta `$ and $`\beta `$, we find these marginal stability curves may vary significantly in reference to those obtained by Shu et al. (2000). Numerical examples are illustrated in Figs. $`510`$ (see also Fig $`D1`$ in Appendix D for a typical solution structure). Generally speaking, the introduction of a gaseous SID tends to decrease the stability of the composite SID system. Meanwhile, the possibility for the system to collapse may be suppressed (Fig. 9). In primordial disk galaxies where gas disk mass exceeds stellar disk mass, the possibility of large-scale collapse may be suppressed while the chance of ring fragmentation increases.

The $`|m|=1`$ case is again similar to the single SID problem. Specifically, the $`y_1`$ solution branch, which is the novel solution of a two-SID system, remains always negative and is thus unphysical. Meanwhile the $`y_2`$ solution branch, which is the counterpart of one-SID problem, can still have physically valid portions constrained by inequalities (102) and (104).

For $`|m|2`$ cases, we have obtained analytical results almost in the same forms of aligned cases. This can be clearly seen by comparing equations $`(94)(99)`$ with equations $`(54)(59)`$. By replacing the Kalnajs function $`𝒩_m(\alpha )`$ with $`(m^2+\alpha ^2+1/4)^{1/2}`$ and $`(m^2+\alpha ^2+1/4)^{1/2}`$ with $`|m|`$ in results of spiral cases, the resulting equations become the same as the aligned cases. To reiterate, the aligned case is just a special one complementary to spiral cases, all of which involve hydrodynamic density waves in a composite system of two coupled SIDs.

It should be mentioned that in the WKBJ regime, such stationary logarithmic spiral density wave perturbations actually carry an outgoing angular momentum flux which is constant in $`r`$. With our composite SID model, a source of angular momentum is thus required at the origin $`r=0`$ where the surface mass density diverges as $`r^1`$. More generally, the outward advective angular momentum flux (stellar and gas SIDs together) is constant in $`r`$, while the gravity torque term involves a nontrivial integral in $`z`$ (see Section 4 of Fan & Lou 1999). In various contexts, several authors, e.g. Lynden-Bell & Kalnajs (1972), Goldreich & Tremaine (1978), and Fan & Lou (1999), have previously investigated angular momentum fluxes carried by density waves or MHD density waves in the tight-winding or WKBJ regime.
3. A Composite System of Two Partial SIDs

In contexts of disk galaxies, the gravitational effect of the more massive dark-matter halo may be incorporated using a model of so-called partial SID by introducing a $``$ parameter which is the ratio of disk potential to the total potential ($`0<1`$). The background surface mass densities of the two coupled SIDs are modified. For stationary perturbations in a composite partial SID system ($`0<<1`$), it is straightforward to make use of the full SID formalism ($`=1`$) by replacing $`1/(1+\delta )`$ with $`/(1+\delta )`$ and $`\delta /(1+\delta )`$ with $`\delta /(1+\delta )`$ in equations (39) and (40).

Naturally, all the previous full SID results are modified by the $``$ parameter quantitatively. In particular, the aligned $`|m|=1`$ case is no longer a solution for arbitrary positive $`D_s^2`$, that is, such stationary eccentric displacements are forbidden in a composite system of two partial SIDs. Also, the critical value $`\beta _c`$ for physically valid $`y_2`$ solution branch (both aligned or unaligned) with $`|m|2`$ is increased somewhat.

While our model formulation is highly idealized in several aspects, it is one step closer to a realistic disk galaxy. The theoretical results here provide a few useful concepts and possible clues for diagnostics of disk galaxies composed of stars and gas coupled by mutual gravity on large scales. Analytical solutions derived here offer valuable benchmarks for numerical codes designed for N-body computations with fluid gas dynamics. The analysis here also paves the way to study theoretical model problems of stationary configurations in a composite system of a stellar SID gravitationally coupled with a gaseous MSID (Lou 2002). It is of considerable interest to extend the current analysis into the nonlinear realm (Galli et al. 2001) to explore the large-scale dynamics of spiral shocks.
4. Stability of a Composite SID System

For both aligned and spiral instabilities, a composite SID system is more complicated than a single SID system in several aspects. Firstly, roughly speaking, the number of possible stationary solutions is doubled except for unphysical $`y_2`$ solutions that must be excluded. Secondly, there are two more dimensionless parameters $`\beta `$ and $`\delta `$ that may cause variations in solution structures. Thirdly, because a composite SID system tends to be less stable, $`D_s^2`$ solutions of lower $`y_2`$ branches may correspond to $`𝒯/|𝒲|`$ ratios considerably less than $`0.14`$ (Ostriker & Peebles 1973). Fourthly, in certain parameter regimes of $`\delta `$ and $`\beta `$ with $`m=0`$, the collapse regime of spiral cases (Shu et al. 2000) may disappear. We shall further examine in detail axisymmetric stability properties of a composite SID system in a separate paper. Fifthly, because of variations in solution structures, one needs to be very careful to determine the most vulnerable configurations with smallest $`D_s^2`$ values. We propose to test the Ostriker-Peebles criterion for a composite disk system by numerical simulation experiments and expect considerably lower $`𝒯/|𝒲|`$ ratio than $`0.14`$ in certain parameter regimes studied here.
5. Dynamics of Magnetized SIDs

For the evolution of an individual spiral galaxy, star formation is one major process (e.g., Kennicutt 1989) and the global star formation rate is a key parameter. In the density wave scenario (Bertin & Lin 1996), spiral arms of higher gas concentration are locations of active star forming regions on smaller scales. Therefore, a protodisk galaxy begins with a gas disk will eventually evolve into a stellar disk. In between the two extremes, the disk galaxy involves a gas disk and a stellar disk with a variable mass ratio. As time goes on, their mass ratio as well as the ratio of velocity dispersions change. At various stages, large-scale instabilities may occur in a composite disk system to alter the evolution track. Furthermore, the presence of magnetic field in the gas disk may give rise to magnetohydrodynamic (MHD) density waves (Fan & Lou 1996; Lou & Fan 1998) and instabilities that can significantly affect the star formation rate.

Magnetic field effects are not included at present merely for simplicity. In disk galaxies, the gaseous disk is magnetized with energy densities of magnetic field, thermal gas and cosmic rays being comparable (Lou & Fan 2003). We derive such information based on synchrotron radio emissions from spiral galaxies. Large-scale magnetic fields typically lie in the disk plane, but there are cases with out-of-plane magnetic fields. On various scales, MHD will play an important role in the gaseous disk and provide valuable dynamical as well as diagnostic information (Fan & Lou 1996; Lou & Fan 1998a; Lou et al. 2001, 2002; Lou 2002). In the current contexts, we simply point out two possible extensions of this theoretical analysis. One is to consider stationary perturbation configurations in a composite partial SID system in which the gas disk is isopedically magnetized (Shu et al. 2000). It is then possible to study the influence of magnetic field on the stability of the system. Such a problem is potentially important to the development of galactic wind models. The other is to consider stationary perturbation configurations in a composite partial SID system in which the gas disk is magnetized with coplanar magnetic fields (Lou 2002; Lou & Fan 2002). The presence of magnetic fields allows for more dynamic freedom in the gas disk and gives rise to fast and slow MHD density waves (Fan & Lou 1996; Lou et al. 2001). We shall pursue these more realistic yet challenging MHD model problems in separate papers.

## Acknowledgments

This research has been supported in part by the ASCI Center for Astrophysical Thermonuclear Flashes at the University of Chicago under Department of Energy contract B341495, by the Special Funds for Major State Basic Science Research Projects of China, by the Tsinghua Center for Astrophysics, by the Collaborative Research Fund from the NSF of China (NSFC) for Young Outstanding Overseas Chinese Scholars (NSFC 10028306) at the National Astronomical Observatory, Chinese Academy of Sciences, and by the Yangtze Endowment from the Ministry of Education through the Tsinghua University. Affiliated institutions of Y.Q.L. share the contribution.

## Appendix A

For the convenience of reference, we here outline a general proof that equation (51) for stationary aligned perturbations always have two branches of real solutions of $`yD_s^2`$ when $`|m|2`$. Note that $`m=0`$ is a rescaling of the background and $`m=1`$ leads to unconstrained $`D_s^2`$. By multiplying a factor $`(1+\delta )`$, the coefficients of quadratic equation (51) can be rewritten in the following forms:

$$C_2^{}=\beta (m^22)(m1)(m+2)(1+\delta ),$$
(124)
$$\begin{array}{cc}\hfill C_1^{}=& 2(m1)[(m+2)\delta (m^22)]\beta \hfill \\ & 2(m^21)[(m^22)\delta +m^2+m2],\hfill \end{array}$$
(125)
$$\begin{array}{cc}\hfill C_0^{}=m(m& 1)[m(m+2)\delta +m^22]\beta \hfill \\ & +2m(m^21)(m\delta +m1),\hfill \end{array}$$
(126)
where the common factor $`(m1)`$ is made explicit in all three coefficients. The determinant $`\mathrm{\Delta }`$ of equation (51) is then given by

$`\mathrm{\Delta }`$ $``$ $`C_1^{}_{}{}^{}24C_2^{}C_0^{}`$
$`=`$ $`4(m^21)^2(𝒜\beta ^2+\beta +𝒞),`$ (127)
where three coefficients $`𝒜`$, $``$, and $`𝒞`$ are defined by

$$𝒜[m^22+(m^2+m2)\delta ]^2,$$
(128)
$$\begin{array}{cc}\hfill =& 2(m^22)(m^2+m2)(\delta ^2+1)\hfill \\ & 2(2m^4+2m^39m^24m+8)\delta ,\hfill \end{array}$$
(129)
$$𝒞=[m^2+m2+(m^22)\delta ]^2.$$
(130)
We can now rewrite equation (A4) in the form of

$$\mathrm{\Delta }=4(m^21)^2\left[𝒜\left(\beta +\frac{}{2𝒜}\right)^2+𝒞\frac{^2}{4𝒜}\right].$$
(131)
Since both $`𝒜>0`$ and

$$𝒞\frac{^2}{4𝒜}=\frac{4m^2(m^22)(m^2+m2)\delta (1+\delta )^2}{[m^22+(m^2+m2)\delta ]^2}>0$$
(132)
for $`m2`$, the determinant of equation (51) is clearly always positive. Therefore, equation (51) will give two branches of real solutions for $`D_s^2=y`$ for $`m2`$.

In the same spirit, we can show numerically that for unaligned or spiral cases of $`|m|1`$, equation (88) always has two branches of real solutions for $`D_s^2=y`$. To avoid repetition and redundancy, we omit that demonstration here.

## Appendix B

We demonstrate here that when $`|m|2`$, the two branches of solution of equations (51) and (88) monotonically decrease with increasing $`\beta `$ as shown in Fig. 1.

Let us consider equation (51) for aligned cases with fixed $`|m|`$ to simplify the computation. For example, with $`|m|=2`$ in equation (51), the three coefficients of equation (51) as modified in Appendix A become

$$C_2^{\prime \prime }=2(1+\delta )\beta ,$$
(133)
$$C_1^{\prime \prime }=(2\delta 1)\beta 3(\delta +2),$$
(134)
$$C_0^{\prime \prime }=(4\delta +1)\beta +3(2\delta +1)$$
(135)
by further removing a common factor 4. Explicitly, the two solutions of equation (51) can then be written as

$$y_{1,2}=\frac{C_1^{\prime \prime }\pm (C_1^{{}_{}{}^{\prime \prime }2}4C_2^{\prime \prime }C_0^{\prime \prime })^{1/2}}{2C_2^{\prime \prime }}.$$
(136)
Substituting expressions $`(B1)(B3)`$ for $`C_2^{\prime \prime }`$, $`C_1^{\prime \prime }`$ and $`C_0^{\prime \prime }`$ into equation $`(B4)`$ and taking the first partial derivative with respect to $`\beta `$, we find that for any positive $`\delta `$, $`y_1/\beta `$ and $`y_2/\beta `$ are always negative. Therefore, both $`y_1`$ and $`y_2`$ solution branches decrease with increasing $`\beta `$. The $`|m|3`$ cases can be proven in a similar manner.

For spiral cases, we follow the same procedure to find that for $`|m|2`$ cases both $`y_1`$ and $`y_2`$ solution branches decrease with increasing $`\beta `$ while for the special case of $`|m|=1`$, the $`y_1`$ and $`y_2`$ solution branches increase and decrease with increasing $`\beta `$, respectively.

## Appendix C

One can readily show that $`\mu ^g/\mu ^s`$ decreases with increasing $`D_s^2`$ for both aligned and unaligned cases when $`|m|2`$ and $`|m|1`$, respectively. For aligned cases, we simply rearrange equation (53) into the following form of

$$\frac{\mu ^g}{\mu ^s}=1\frac{(m^22)(1+\delta )}{|m|}+\frac{2(m^21)(1+\delta )}{|m|(D_s^2+1)}$$
(137)
which shows obviously that $`\mu ^g/\mu ^s`$ decreases with increasing $`D_s^2`$ when $`|m|2`$.

For spiral cases, we can also arrange equation (92) into the parallel form of

$$\begin{array}{cc}\hfill \frac{\mu ^g}{\mu ^s}=& 1\frac{(m^22)(1+\delta )}{𝒩_m(\alpha )(m^2+\alpha ^2+1/4)}\hfill \\ & +\frac{[m^22+(m^2+\alpha ^2+1/4)](1+\delta )}{𝒩_m(\alpha )(m^2+\alpha ^2+1/4)(D_s^2+1)}\hfill \end{array}$$
(138)
which again shows clearly that $`\mu ^g/\mu ^s`$ decreases with increasing $`D_s^2`$ when $`|m|1`$ for any radial wavenumber $`\alpha `$. Note in particular that when $`|m|=1`$, the coefficient of $`1/(D_s^2+1)`$ on the right-hand side of equation $`(C2)`$ is positive.

## Appendix D

In general, for two similar dynamical systems coupled together, we expect two sets of oscillation modes (e.g., Vandervoort 1991a, b). This should be true for the current problem of a composite SID system as compared to a single SID system. For the cases of marginal stability for spiral configurations with $`m=0`$, we thought that there might be something interesting from the additional solution of a composite SID system in contrast to the single SID case (Shu et al. 2000). It turns out that this additional solution of $`D_s^2`$ is negative by numerical exploration. The numerical example of Figure $`D1`$ shows the typical structure of solution curves for stationary logarithmic spiral configurations when $`m=0`$, $`\delta =0.2`$ and $`\beta =1.5`$. For mathematical completeness, the negative solution branches are shown explicitly. In Fig. 5 of the main text, the negative solutions of $`D_s^2`$ are discarded.