Global Perturbation Configurations in a Composite Disc System with an Isopedic Magnetic Field

We construct stationary global configurations of both aligned and unaligned logarithmic spiral perturbations in a composite disc system of stellar and isopedically magnetized gaseous singular isothermal discs (SIDs) coupled by gravity. The thin gaseous SID is threaded across by a vertical magnetic field B_z with a constant ratio of the surface gas mass density to B_z. There exist two classes of stationary magnetohydrodynamic (MHD) solutions with in-phase and out-of-phase density perturbations. For both aligned and unaligned cases with azimuthal periodicities |m|\geq 2 (m is an integer), there may be two, one, and no solution situations, depending on the chosen parameters. For the transition criteria from an axisymmetric equilibrium to aligned secular bar-like instabilities, the corresponding T/|W| ratio can be much lower than the oft-quoted value of T/|W|\sim 0.14, where T is the total rotational kinetic energy and W is the total gravitational potential energy plus the magnetic energy. The T/|W| ratios for the two sets of solutions in different ranges are separated by m/(4m+4). For the unaligned cases, we study marginal stabilities for axisymmetric (m=0) and non-axisymmetric (m\neq 0) disturbances. The gravitational influence of an axisymmetric dark matter halo on the background is also examined. The global analytical solutions and their properties are valuable for testing and benchmarking numerical MHD codes. Our model analysis contains more realistic elements and offers useful insights into the structures and dynamics of disc galaxies consisting of stars and magnetized interstellar medium. In the presence of star burst activities, supernovae, hypernovae, superbubbles etc., our open magnetic field geometry in disc galaxies bears strong implications on circumnuclear and spiral galactic winds.


INTRODUCTION
The theoretical magnetohydrodynamic (MHD) disc model we set out to formulate in this paper is to explore possible large-scale global perturbation structures and stationary MHD density waves (Fan & Lou 1996;Lou & Fan 1998) in a composite system of a stellar disc and an isopedicallly magnetized gaseous disc intended for the interstellar medium (ISM). These two gravitationally coupled discs are approximately treated as 'fluid' and 'magnetofluid' respectively, and are both geometrically idealized as razor-thin singular isothermal discs (SIDs) with the gaseous SID being threaded across by an almost vertical magnetic field throughout. In astrophysical contexts of large-scale structures in disc galaxies, we also include gravitational effects of a massive axisymmetric dark matter halo and adopt a background composite system of two coupled partial SIDs (Syer & Tremaine 1996;Shu et al. 2000;, 2004aLou & Zou 2004Shen, Liu & Lou 2005). Our motivation is to construct solutions with combined analytical and numerical techniques, to understand their basic properties, to provide observational diagnostics, and to reveal or speculate physical implications.
Chakrabarty, Laughlin & Shu (2003) studied substructures in grand-design spiral galaxies, such as branches, spurs and feathers, using a two-component disc model in which the gas component responds passively and nonlinearly to the potential of a rigidly rotating spiral structure involving old stars and halos. We here treat a dynamically coupled twocomponent disc system without or with a massive dark matter halo and focus on stationary global MHD density wave configurations. Specifically, we construct stationary MHD configurations for aligned and unaligned logarithmic spiral perturbations in a composite disc system of two SIDs with flat rotation curves. For observational diagnostics of nearby disc galaxies, we also examine phase relationships among perturbation patterns of the stellar surface mass density, the ISM surface mass density and the isopedic magnetic field.
We now proceed to provide the more general background information relevant to the idealized MHD composite disc problem to be formulated and investigated here.
In a pioneering work on a composite disc system of stellar and gaseous discs dynamically coupled by gravity, Lin & Shu (1966, 1968) proposed a combined approach involving a distribution function for the stellar disc and a fluid description for the gas disc to derive and analyze the local dispersion relation of coplanar galactic spiral density waves in the WKBJ approximation. The basic physical scenario is that stars form out of gas clouds in the ISM disc, leading to the coexistence of a stellar disc and a magnetized ISM disc at a later evolution stage of disc galaxy. Since the seminal work of Lin & Shu, there have been extensive theoretical researches on density wave oscillations, perturbation configurations and stability properties of a rotating composite disc system, mainly in the galactic context. Kato (1972) studied oscillations and overstabilites of density waves using a formalism similar to that of Lin & Shu (1966, 1968. Jog & Solomon (1984a, b) discussed the growth of local axisymmetric perturbations using a two-fluid formalism in a composite disc system. Bertin & Romeo (1988) investigated the role of a gas disc for spiral modes in a two-fluid model framework. The influence of interstellar gas on oscillations and stabilities of spheroidal galaxies was studied by Vandervoort (1991a, b). In order to account for the effects of the disc thickness, Romeo (1992) adopted a two-fluid approach to investigate a two-component disc system with finite disc thickness. Lowe et al. (1994) performed an extensive analysis for morphologies of disc galaxies. Different effective Q ef f parameters (Safronov 1960;Toomre 1964) have been suggested for the axisymmetric stability of a composite disc system in a two-fluid formalism by Elmegreen (1995) and Jog (1996). Lou & Fan (1998b) recently used the two-fluid formalism to study properties of open and tight-winding spi-ral density-wave modes in a composite disc system.  discussed stationary global perturbation structures in a two-fluid formalism and offered a more straightforward D−criterion for the axisymmetric instability in a composite SID system instead of a redefinition of a new Q ef f parameter . Considering the magnetic field embedded in the gaseous disc, Lou & Zou (2004 studied the stationary global coplanar MHD perturbation structures and axisymmetric stability in a composite SID system.
A series of astrophysical disc problems involves stability analysis of a SID. Since the pioneer work of Mestel (1963) more than four decades ago, numerous theoretical and numerical studies have been carried out in this subject area (e.g. Zang 1976;Toomre 1977;Lemos, Kalnajs & Lynden-Bell 1991;Lynden-Bell & Lemos 1999;Goodman & Evans 1999;Chakrabarti, Laughlin & Shu 2003). An important breakthrough by Syer & Tremaine (1996) was to derive the semi-analytic solutions for stationary coplanar perturbation configurations in a class of power-law discs (i.e., SID is only a special case). Shu et al. (2000) obtained stationary solutions for global perturbation configurations in an isopedically magnetized SID with a flat rotation curve. Based on extensive numerical calculations, they interpreted these stationary aligned and unaligned logarithmic spiral configurations as the onsets of bar-type and barred-spiral instabilities (see also Galli et al. 2001). As a different yet complementary work to the analysis of Shu et al. (2000),  studied a coplanar MHD perturbation analysis in a single background SID embedded with an azimuthal magnetic field, from the perspective of stationary fast and slow MHD density waves (FMDWs and SMDWs; Fan & Lou 1996;Lou & Fan 1998a).  also derived a form of magnetic virial theorem for an MSID and suggested by analogy that the ratio of the disc rotational kinetic energy to the sum of the gravitational and magnetic energies may be crucial for the onset of bar-like instability in an MSID system. To model large-scale effects of a galactic magnetic field, it would be more realistic to treat large-scale structures and dynamics as a composite disc system of stars and a magnetized ISM.  made an foreray into this model problem, constructed stationary aligned and unaligned logarithmic spiral configurations in such a composite SID system and further examined axisymmetric instability properties . Meanwhile, Lou & Zou (2004 incorporated a coplanar magnetic field into the composite SID system for a more specific analysis. The geometry of an isopedic magnetic field proposed by Shu & Li (1997) is an open magnetic field configuration across the disc plane. For disc galaxies, open magnetic fields are most likely interlaced or intermingled with 'closed' coplanar magnetic fields in a stochastic manner. ⋆ For a much simplified theoretical model analysis at this stage, we take an isopedic magnetic field geometry all over a composite SID system. Together with other physical mechanisms (e.g., massive star formation activities, supernovae, hypernovae, hot bubbles and superbubbles etc.), such an open magnetic field geometry can generate and channel hot gas outflows, often referred to as 'galactic winds' under increasing scrutiny of multi-band observations and the X-ray band in particular. Along with observational and theoretical studies of stellar winds (e.g., Burke 1968;Holzer & Axford 1970), an initial impetus was made to formulate the seminal concept of galactic-scale winds blowing away from both sides of the disc plane (e.g., Johnson & Axford 1971). As galactic winds are expected to be a relatively weak phenomenon in terms of radiative signatures, we start to accumulate some substantial observational evidence in recent years. It is generally conceived that supernova explosions and stellar radiations drive galactic outflows. Strickland, Ponman & Stevens (1997) have reported X-ray observations of galactic outflows from the starburst galaxy M82, which may imply an overall configuration of a galactic wind. On both observational and theoretical grounds, magnetic field configurations away from a galactic plane should play dynamical, diagnostic and geometric roles in probing galactic winds. As a result of such open magnetic field, the outflow of hot gas can lead to largescale observable effects.
For such a magnetized composite disc system, MHD will play an indispensable role in the gaseous MSID and reveal more realistic aspects of coupled large-scale dynamics as well as diagnostic information. These important MHD disc problems (Shu et al. 2000;Lou & Zou 2004Shen, Liu & Lou 2005) are not only interesting for their own sakes, but also serve as necessary steps for developing more and more realistic physical models. We construct here stationary perturbation configurations for aligned and unaligned logarithmic spiral cases in a composite system of a stellar SID and an isopedically magnetized gaseous SID, and discuss their stability properties. For numerical MHD simulations into the nonlinear regime, the global perturbation solutions here are valuable for initializing, benchmarking and understanding the MHD code development.
We adopt a relatively simple formalism for a composite disc system containing infinitely thin fluid (for stars) and magnetofluid (for the magnetized ISM) discs dynamically coupled by gravity. We describe in Section 2 the isopedic magnetic field and a few important theorems for the gaseous MSID. In Section 3, we determine conditions for the background axisymmetric equilibrium state for both stellar SID and gaseous MSID, and derive linearized equations for MHD perturbations. The global MHD perturbation solutions can be classified as aligned perturbations and unaligned logarithmic spiral perturbations; these solutions are analyzed in detail in Section 4. The exact solutions of stationary MHD perturbations, their properties and their corresponding phase relationships among perturbation variables are examined and summarized in Section 5. Several details of mathematical analysis are included in Appendices A and B for the convenient of references.

ISOPEDIC MAGNETIC FIELDS
Magnetic fields may be generated and sustained through nonlinear MHD dynamo processes in an electrically conducting gas medium in motion (Parker 1979 and extensive references therein). Several numerical simulation results and some theories for galaxy formation lend supports to the basic assumption that the mass-to-magnetic flux ratio Λ may be constant in space and time. For example, model calculations and numerical simulations of cloud core formation by ambipolar diffusion tend to produce gravitationally unstable inner regions with approximately constant values of Λ (e.g., Nakano 1979;Lizano & Shu 1989). In the models of Basu & Mouschovias (1994), the mass-to-magnetic flux ratio Λ remains almost constant when the variation ranges of density and magnetic field flux are large. In the theoretical formulation of Shu & Li (1997) and Shu et al. (2000), it was thus presumed that Λ is a constant; in particular, they referred to these earlier results to justify this assumption. Here, we emphasize the theoretical fact that Λ remains constant is a natural consequence of the standard ideal MHD equations. In cylindrical coordinates (r, ϕ, z) and by using the gas mass conservation equation (20) below and theẑ−component of the magnetic induction equation in the thin disc geometry, one can readily show that where the relevant notations bear the conventional meanings in MHD. In reference to the works of Mouschovias (1994) and Li & Shu (1996) as well as the earlier numerical results (Nakano 1979;Lizano & Shu 1989), it is fairly clear that the above equation provides a natural explanation for the fact that Λ remains constant in the evolution of an MHD disc. For the model problem under consideration, we ignore vertical structures within the disc thickness as the horizontal scale of the problem is very much larger. To analyze properties of a drastically flattened disc, we adopt the cylindrical coordinates (r, ϕ, z) in our mathematical model prescription. We assume that the mass density ρ and the electric current density je in our model to have nonzero values only within a narrow range ∆z ≪ r about the disc midplane z = 0. We are justified to define a surface mass density Σ and a surface electric current density J by Σ(r, ϕ, t) ≡ ∆z ρ(r, ϕ, z, t)dz and J(r, ϕ, t) ≡ ∆z je(r, ϕ, z, t)dz , with J = Jrêr + Jϕêϕ having only two orthogonal components coplanar with the disc and the z−component of J within the disc thickness is omitted.
As for the two SIDs, one is a stellar SID approximated as a fluid for large-scale dynamics and the other is a magnetized gas SID treated as a magnetofluid. Over large scales, we ignore the interaction between the magnetic field em-bedded in the ISM and the stellar disc. It then comes the physical scenario that the stellar and magnetized gas discs are dynamically coupled by gravity, while the magnetic field directly interacts with the gaseous ISM disc.
On the basis of the Poisson equation relating the gravitational potential and the mass density, there will also be a jump in the vertical gravitational field given by 4πGΣ across a thin disc. In other words, the vertical gravitational field just above and below the disc point towards the midplane z = 0 and have the same strength of 2πGΣ. We choose to define the ratio η of the horizontal gravitational acceleration f (continuous across the thin disc in vertical direction) to the vertical gravitational acceleration just above the disc as a dimensionless parameter, namely and the local value of the mass-to-magnetic flux ratio as We then take Λ to be constant in both space and time, when computed at the footpoint of a magnetic field line anchored through the disc (the adjective 'isopedic' comes from the Latin word 'ped' for foot). For ideal MHD with the above frozen-in condition, an initially isopedic disc retains the same value of Λ throughout the entire MHD evolution. It is convenient and advantageous to nondimensionalize this ratio in the form of We now discuss the two MHD theorems proven by Shu & Li (1997) that form the fundamental part of the MHD equations to come in our theoretical model analysis.
THEOREM 1: For an arbitrary distribution of the surface mass density Σ(r, ϕ, t), which may or may not be in mechanical equilibrium in the horizontal directions, the magnetic tension force in the plane of a thin isopedic disc equals −1/λ 2 times the horizontal self-gravitational force f .
THEOREM 2: If we also assume that the disc remains in a vertical magnetostatic equilibrium, then the magnetic pressure integrated over the disc thickness may be approximated as (1 + η 2 )/(λ 2 + η 2 ) times the thermal gas pressure integrated over z.
By the first theorem, the magnetic tension force ften tangential to the disc plane is simply ften = −f /λ 2 and the combined forces of magnetic tension and gravitational field act in a form of a diluted horizontal gravity where ǫ ≡ 1 − 1/λ 2 is the reduction factor.
By the second theorem, it follows approximately that where Πm is the projected magnetic pressure acting on the gas material and Πg is the vertically integrated thermal gas pressure. Therefore, the total projected pressure is where the enhancement factor Θ is given by Coming to a composite disc system of two gravitationally coupled discs with the gas SID being isopedically magnetized, the above two theorems should be modified accordingly. By carefully examining the two theorems in details, we note that as the stellar disc does not interact with the magnetic field, the first theorem of Shu & Li should retain its basic form with some relevant parameters being properly adjusted. For example, Λ should now stand for Σ g /Bz where Σ g is the surface mass density of the gas disc and the magnetic tension force ften acts on the magnetized gas SID. Naturally, we introduce an appropriate η parameter for convenience, namely η ≡ f g /(2πGΣ g ), where f g is the contribution to the total f from the gaseous MSID. For the two SIDs having the same structure and density profile, we further obtain where Σ s is the surface mass density of the stellar SID, Σ is the total surface mass density and f s is the contribution to the total f from the stellar SID. For the second theorem of Shu & Li, the generalization and extension are somewhat different. We now examine the procedure of Shu & Li to generalize the second theorem for a composite disc system with one stellar SID and one isopedically magnetized SID.
The force per unit volume associated with the magnetic pressure gradient is given by By integrating −∇[B 2 /(8π)] from z = −∞ to z = +∞, the contribution from the z−component of the magnetic pressure gradient vanishes for B 2 → 0 as |z| → +∞. We thus obtain the parallel magnetic force per unit projected surface area as where the parallel gradient operation within the disc plane is defined by Away from the disc, the total integral of equation (3) does not act on materials contained in the gaseous MSID. In the vacuum regions above and below the MSID, the horizontal force of magnetic pressure gradient does not vanish but counteracts instead a nonzero force of magnetic tension, such that the sum produces a force-free environment in vacuum (i.e., no electric current flowing around). We specifically denote the part of the projected magnetic pressure acting on the gaesous MSID as Πm. For a magnetized gas disc of a thin but nonvanishing thickness, we estimate Πm by where z0 is the effective half-thickness of the gaseous MSID and B 2 + is the value of B 2 evaluated just above the upper disc surface. While our estimate for Πm is plausible on the dimensional ground, the possibility remains for a multiplicative factor of order unity whose exact value would depend on the vertical structure within a thin MSID.
In order to determine the two expressions of z0 and B 2 + , we begin with a generalization which states that a quasi vertical magnetostatic equilibrium in a thin gas disc requires the sum of the gas, magnetic and gravitational pressures to be independent of the vertical height z, namely, where Pg is the thermal pressure of the gaseous MSID and K is an integration constant; and by assuming that the gaseous MSID and the stellar SID interact only via the mutual gravity, we have the total effective surface mass density in equilibrium as In the midplane z = 0, we have by symmetry σ = 0, Br = Bϕ = 0, Bz = Σ g /Λ, while Pg = Πg/(2z0) with Πg being the vertically integrated thermal gas pressure and z0 is the half-scaleheight in the vertical direction. Thus equation (5) yields at z = 0 the following expression for constant K, Meanwhile for r ≫ z ≫ z0, we have Pg = 0, Ps = 0, σ = Σ s + Σ g , and B = (Σ g /Λ)êz − f g /(2πGΛ) [see derivations of Shu & Li (1997) that lead to their equations (2.14) and (2.15)]. In these spatial regions, equation (5) leads to By introducing the gas-to-stellar surface mass density ratio δ ≡ Σ g /Σ s from the two SIDs, we derive from equations (7) and (8) the following relation It then follows immediately that Just above the disc, equation (5) leads to Using expression (7) for constant K, we then obtain Multiplying equation (12) by 2z0 with z0 being given explicitly by expression (10), we derive from definition (4) It is convenient to introduce a modified λ parameterλ ≡ (1 + δ −1 )λ such that expression (13) becomes By summing up the gas and magnetic pressure contributions, we obtain the total pressure in the form of where the modified enhancement factor Θ is given by This completes our extension and generalization of the second theorem of Shu & Li as applied to the MHD problem of a composite disc system consisting of one stellar SID and one gaseous isopedic MSID. At the end of this section, we summarize the two extended theorems for a composite MSID system with an isopedic magnetic field as follows. In the presence of an isopedic magnetic field, the horizontal self-gravity force f g and the gas pressure integral Πg(= a 2 g Σ g ) of the gaseous MSID will be modified as where ǫ ≡ 1 − λ −2 and Θ is defined above by equation (16) that involves essentially three dimensionless parameters δ ≡ Σ g /Σ s , λ ≡ 2πG 1/2 Σ g /Bz, η ≡ f g /(2πGΣ g ) = f /(2πGΣ) andλ ≡ (1 + δ −1 )λ.

Basic Equations for a Composite MSID System and the Background Equilibrium
The two coupled SIDs located at z = 0 are both approximated as infinitesimally thin for expediency and they are coupled through the mutual gravity. For large-scale stationary MHD perturbations, diffusive processes such as viscosity, ambipolar diffusion and thermal diffusion etc. are all ignored in our formulation for simplicity. For physical variables under consideration, we use either superscript or subscript s to indicate the association with the stellar SID and either superscript or subscript g to indicate the association with the gaseous MSID. In the cylindrical coordinates (r, ϕ, z), the ideal fully nonlinear fluid-magnetofluid equations for a composite MSID system can be readily written out. For the stellar SID in the 'fluid' approximation, we have for the two-dimensional stellar mass conservation, for the azimuthal momentum equation, where the 'isothermal' approximation has been invoked. In parallel, for the gaseous isopedic MSID in the magnetofluid approximation, we have for the two-dimensional gas mass conservation, for the radial momentum equation, and ∂j g ∂t for the azimuthal momentum equation, where ǫ ≡ 1 − 1/λ 2 and the modified Θ here is defined by equation (16). The dynamical coupling between the two sets of fluid and magnetofluid equations is caused by the gravitational potential through the Poisson integral relations and as well as through the effects of δ and η in the Πm term. In equations (17)−(24), Σ s is the stellar surface mass density, u s is the radial component of the bulk stellar 'fluid' velocity, j s is the z−component of the specific angular momentum of the stellar 'fluid', and as is the stellar velocity dispersion presumed to be constant (or an equivalent 'isothermal sound speed' for the stellar SID), a 2 s Σ s stands for an effective pressure in the isothermal approximation. For physical variables of the gaseous isopedic MSID, we simply replace the superscript s by g systematically. Here, we assume that the stellar and magnetized gaseous SIDs interact mainly through the mutual gravitational coupling on large scales. As the isopedic magnetic field only interacts with the gas disc, Θ and ǫ only appear in the MHD equations for the gaseous MSID.
Before an MHD perturbation analysis, one needs to establish a background magneto-rotational equilibrium for the composite MSID system that is dynamically consistent with equations (17)−(24). For that purpose, we presume an axisymmetric background for both stellar and magnetized gaseous SIDs, with the same form of power-law surface mass densities (i.e., Σ ∝ r −1 ) yet with different flat rotation curves in general.
Using the radial momentum equations (18) and (21) for the background equilibrium of axisymmetry with u s 0 = u g 0 = 0, Ωs = j s 0 /r 2 , and Ωg = j g 0 /r 2 , we readily obtain and −Ω 2 where δ ≡ Σ g 0 /Σ s 0 , and Ωs ≡ asDs/r and Ωg ≡ Θ 1 2 agDg/r for the two rotational Mach numbers Ds and Dg, respectively. The two equilibrium surface mass densities Σ s 0 and Σ g 0 can be immediately determined as and For stationary global MHD perturbation configurations to persist, Ds and Dg can only take on specific values (e.g. , 2004aLou & Zou 2004. In the above equations, the two parameters δ and ǫ are constant, while Θ is a variable in general. However, in our MHD perturbation analysis, Θ variation will be ignored merely for simplicity without losing essential physics (see Shu & Li 1997;Shu et al. 2000). We now further introduce a new dimensionless parameter β ≡ a 2 s /a 2 g for the square of the 'sound speed ratio'.
Because |f g | = |dφ g 0 /dr| = 2πGΣ g 0 in an axisymmetric equilibrium state, we take η = 1 and have Θ > 1 in general. For a positive Σ g 0 as a necessary physical requirement, the equilibrium solution (28) implies the following three inequalities in sequence As δ > 0 by definition, it is obvious that λ(δ −1 + 1) > λ. A multiplication of this inequality with the third inequality above on both sides leads toλ ≡ λ(1 + δ −1 ) > 1 for a physical magneto-rotational background equilibrium. It then follows that the enhancement factor Θ defined by equation (16) should be constrained within a finite range of 1 < Θ < 2. We now discuss the parameter ǫ. In Shu & Li (1997), they require ǫ ≥ 0 which is the same as λ ≥ 1 for a gravitationally bound disc. In the case of a single SID, the physical explanation of this condition is that a rotating disc with self-gravity must be held together in a radial force balance. As a result of rotation and singular isothermal density distribution, the centrifugal force and the effective and thermal pressure forces all point radially outward. The only counterforce to keep the disc in equilibrium is the radially inward gravity force. For ǫ < 0, the magnetic tension force will be stronger than the horizontal gravitational force; in this case, there will be no rotational equilibrium for the disc system at all to begin with. In our composite model containing one SID and one isopedic MSID, the necessary condition for an equilibrium to exist is 1 + ǫδ > 0 in comparison with the condition of ǫ > 0 of Shu & Li. More precisely, we should have 0 < ǫ ≤ 1 for a single SID and 0 < (1 + ǫδ) ≤ (1 + δ) for a composite system of two gravity-coupled SIDs.

MHD Perturbation Equations
We now introduce nonaxisymmetric MHD perturbations to the disc system, denoted by physical variables with subscript 1, in both stellar SID and gaseous isopedic MSID, namely It is then straightforward to write down the linearized MHD perturbation equations in the forms of ∂j s 1 ∂t + rΩsu s 1 + Ωs for the stellar SID, and ∂j g 1 ∂t + rΩgu g 1 + Ωg for the gaseous isopedic MSID, with the total gravitational potential perturbation φ1 = φ s 1 + φ g 1 given by the sum of Assuming a Fourier harmonics in the periodic form of exp[i(ωt − mϕ)] for MHD perturbation solutions in general (after taking the real part), we write for perturbations where l = s or g denotes associations with stellar or gaseous discs respectively. By substituting expressions (37)−(40) into equations (29)−(36), we derive for the stellar SID for the gaseous isopedic MSID and for the two gravitational potential perturbations We use equations (42) and (43) to express U s and J s in terms of Ψ s ≡ (a 2 s µ s /Σ s 0 ) + V s + V g for the stellar SID and similarly, use equations (45) and (46) to express U g and J g in terms of Ψ g ≡ (Θa 2 g µ g /Σ g 0 ) + V s + ǫV g for the gaseous MSID. In the following analysis, we shall assume for simplicity that Θ remains constant in the presence of MHD perturbations in a composite MSID system (see Shu et al. 2000). The resulting expressions then become and for the stellar SID, and and for the gaseous isopedic MSID, respectively. A substitution of expressions (49) and (50) into equation (41) leads to for the stellar SID. Likewise, a substitution of expressions (51) and (52) into equation (44) leads to for the gaseous isopedic MSID. Equations (53) and (54) should be solved together with Poisson integral relations (47) and (48). For stationary MHD perturbation solutions (ω = 0) with zero pattern speed, we readily have for the stellar SID, and for the gaseous isopedic MSID, respectively. Using conditions and relationships (25)−(28), we further reduce the above two equations to the forms of for the stellar SID and the gaseous isopedic MSID, respectively. For nonzero Ωs and Ωg, we have for the stellar SID and for the gaseous isopedic MSID, respectively. Equations (59) and (60) are to be solved simultaneously with Poisson integrals (47) and (48) in order to determine the two rotational Mach numbers D 2 s and D 2 g .

A Composite System of Partial (M)SIDs
Since the early pioneer work of Zwicky (1933) and Smith (1936), evidence has gradually emerged for the presence of dark matters in clusters of galaxies. In the 1970s, systematic observations of field spiral galaxies further revealed the generic presence of dark matters in isolated spiral galaxies.
By careful and independent measurements of galactic rotations (e.g., Rubin 1987;Kent 1986Kent , 1987Kent , 1988, the more or less flat rotation curves of galaxies determined by different methods clearly imply the presence of dark matter halos in spiral galaxies. To make our formulation more general and realistic, we here propose a theoretical model of a composite partial MSID system to take into account of the Newtonian gravitational effect of a presumed axisymmetric dark matter halo (e.g., Binney & Tremaine 1987;Bertin & Lin 1996;Syer & Tremaine 1996;Shu et al. 2000;, 2004aLou & Zou 2004Shen, Liu & Lou 2005). For a composite partial MSID system, the preceding mathematical formulation for a composite MSID system remains to be its essential part yet with an important distinctive feature of adding an axisymmetric gravitational potential in the background magnetorotational equilibrium. In numerical simulations designed for galaxy formation, velocity dispersions of dark matter 'particles' in the halo are typically high (e.g., Hoeft et al. 2004).
For simplicity, we may therefore ignore dynamical perturbation responses of this massive dark matter halo to MHD perturbations in the composite partial MSID system. More specifically, we introduce an additive gravity term −∂Φ/∂r to the basic equations (18), (19) and (21), (22), where Φ represents the gravitational potential contribution from the background dark matter halo of axisymmetry. In the two components of the momentum equation for the stellar disc, we then have meanwhile, for the two components of the momentum equations in the gaseous isopedic MSID, we have For such a composite partial MSID system, we now introduce a dimensionless ratio F parameter defined by F≡ (φ g +φ s )/(φ g +φ s +Φ) Lou & Zou 2004. The background magneto-rotational equilibrium should be modified accordingly. As before, we still write Ωs = asDs/r, Ωg = agDg/r, and u s 0 = u g 0 = 0 for the background. The background equilibrium equations then become It follows immediately that the two modified equilibrium surface mass density profiles are now given by and respectively. From these two radial equilibrium conditions, we readily derive the requirement ) resulting from the same total gravitational potential that couples the two SIDs. We refer to such a system as a composite system of partial (M)SIDs for 0 < F < 1, in contrast to a composite system of full (M)SIDs with F = 1.
The linearized MHD perturbation equations in a composite system of partial (M)SIDs should take the same forms as those in a composite system of full (M)SIDs explicitly written down earlier in this section. For stationary perturbation configurations in the modified equilibrium with ω = 0, we then have for the stellar partial SID and simultaneously for the gaseous partial MSID where F factor appears explicitly in both relations.

Two Different Limiting Procedures
Here, we discuss specifically the axisymmetric case of m = 0 and start from the very beginning with ω = 0 in Fourier harmonics (37)−(40). By substituting expressions (37)−(40) with m = 0 into equations (29)−(36), we derive for the stellar SID, for the gaseous isopedic MSID, and for the two gravitational potential perturbations. Again, we use equations (73) and (74) to express U s and J s in terms of Ψ s ≡ (a 2 s µ s /Σ s 0 ) + V s + V g for the stellar SID and in parallel, we use equations (76) and (77) to express U g and J g in terms of Ψ g ≡ Θ(a 2 g µ g /Σ g 0 ) + V s + ǫV g for the gaseous MSID. The resulting expressions are for the stellar SID, and for the gaseous MSID, respectively. A substitution of expressions (80) and (81) into equation (72) leads to for the stellar SID. Likewise, a substitution of expressions (82) and (83) into equation (75) leads to for the gaseous isopedic MSID. Ordinary differential equations (ODEs) (84) and (85) are to be solved with Poisson integrals (78) and (79) simultaneously. The case of m = 0 is special. If we set ω = 0, all MHD perturbation equations are satisfied in a trivial manner. The harmonic factor exp [i(ωt − mϕ)] becomes unity for both ω = 0 and m = 0. In order to construct stationary configurations, we set ω → 0 to derive limiting solutions.
The two equations above can be cast into 86) for stellar disc and for the gaseous isopedic MSID. For small ω = 0, we readily obtain and by removing ω = 0 factor in equations (86) and (87). For the onset of marginal instability in the limit of ω → 0, we then derive from equations (88) and (89) and These two resulting equations are exactly the same as those obtained by setting m = 0 in equations (57) and (58), for the stellar SID, and for the gaseous isopedic MSID, respectively. This situation of our composite MSID system parallels the case of Shu et al. (2000) for a single SID which is isopedically magnetized as well as the case of  for a composite SID system. For a background MSID with a coplanar magnetic field, the two different limiting procedures of m = 0 and ω → 0 give rise to two different sets of resulting differential equations Lou & Zou 2004;Shen, Liu & Lou 2005).

Aligned Global Perturbation Configurations
By definitions, aligned MHD disturbances in a disc mean that maxima and minima of perturbations occurring at different radii line up in the azimuth. In contrast, when maxima and minima of perturbations in a disc shift systematically in the azimuthal angle from one radial location to the next, the resulting configuration for such disturbances are referred to as unaligned or spiral patterns (Kalnajs 1973;Binney & Tremaine 1987;Shu et al. 2000;Lou & Zou 2004).

Axisymmetric Disturbances (m = 0)
With ω = m = 0, the solution to equations (41)−(48) takes the form of U s = U g = 0, µ s = K s 1 /r, µ g = K g 1 /r, J s = K s 2 r, J g = K g 2 r, V s = K s 3 ln r, and V g = K g 3 ln r, where K s,g i for i = 1, 2, 3 are all constants which can be chosen properly such that equations (42), (45), (47), (48) are satisfied. Nevertheless, such a solution merely represents a rescaling of one axisymmetric equilibrium to a neighbouring axisymmetric equilibrium. Such rescaling solutions are allowed but are not of particular interest here.

Cases with |m| ≥ 1: Nonaxisymmetric Disturbances
In discs with power-law distributions of background physical variables, we consider aligned perturbation variables that carry the same power-law dependence as those of the equilibrium SID do. By this choice, we make use of the following exact potential-density pair relation, namely where σs and σg are two small constants for perturbations. A substitution of equations (94) and (95) into equations (59) and (60) then leads to the following equations: where the four coefficients H1, H2, G1 and G2 are defined by By substituting the forms of µ s and µ g given by equation (94) into equations (96) and (97), we readily obtain (1 − H2m 2 )µ g = G2m 2 µ s .
For nontrivial solutions of µ s and µ g , we then derive for the stationary dispersion relation to determined D 2 s and thus D 2 g . Using definitions (98)−(101), the solution condition (104) can be written explicitly in the form of In terms of background profiles (27) and (28), we obtain and then where the ratio β ≡ a 2 s /a 2 g was introduced earlier. As is obvious, equation (105) remains unchanged by replacing m with −m. For simplicity, we use m to represent |m|. In fact, positive and negative m correspond to trailing and leading spiral patterns relative to the sense of SID rotation, respectively. After straightforward manipulations and rearrangements, equation (105) can be cast into a quadratic equation in terms of y ≡ D 2 s , namely where the three coefficients C2, C1 and C0 are respectively. From these three coefficient expressions, the resulting y solutions remain valid for both positive and negative signs of m. Unless otherwise stated, we shall take m to be positive without loss of generality. As seen from our MHD perturbation analysis so far, the Θ parameter always appears together with the β parameter; it is natural and convenient to introduce a new parameter χ ≡ β/Θ to simplify notations. Physically, Θa 2 g represents a sort of magnetosonic speed squared given the isopedic magnetic field geometry; the ratio χ ≡ a 2 s /(Θa 2 g ) thus stands for the square of the ratio between the stellar velocity dispersion and the magnetosonic speed.
We now discuss properties of quadratic equation (108) with its determinant ∆ given explicitly by Except for the case of m = 1, the important fact that ∆ > 0 for m > 1 means that quadratic equation (108) always have two real solutions of y ≡ D 2 s . The physical meaning is that for aligned nonaxisymmetric stationary configurations to exist in a composite MSID system with an isopedic magnetic field, condition (105), or equivalently, condition (108) must be satisfied for appropriate values of D 2 s given a set of specified parameters m, δ, ǫ and χ ≡ β/Θ. The two real y ≡ D 2 s solutions may not be necessarily physical unless the nonnegative requirements of both D 2 s ≥ 0 and D 2 g ≥ 0 are met. We now elaborate several subtle points below. For the aligned case with m = 1, it turns out that C2 = C1 = C0 = 0 by definitions. Equation (108) can therefore be satisfied for arbitrary values of D 2 s . This situation is quite similar to the m = 1 cases of a single isopedically magnetized SID studied by Shu et al. (2000) and of a composite system of two coupled SIDs analyzed by .
When m ≥ 2, the expressions of coefficient parameters in the equation can be rearranged into the following forms by multiplying through a factor of (1 + δ) 2 , namely Since 1 + δ > 0, 1 + ǫδ > 0 and χ > 0, we further simplify these three coefficient expressions to obtain equivalently For the convenience of further analysis, we introduce two new dimensionless parameters B and X defined by It then follows that C2, C1 and C0 can be expressed as Because there are several parameters in our theoretical model analysis and our main motivation is for galactic applications, it would be efficient and more sensible to have rough ranges to bracket these parameters in order to explore different regimes for various possible solutions numerically.
For this purpose, we estimate ranges of the parameters λ, β, ǫ and Θ for disc galaxies. A typical magnetic field in a spiral galaxy is about 1 ∼ 10µG and higher values of ∼ 30 − 40µG may be reached in central circumnuclear regions by equipartition arguments (e.g., Lou et al. 2001). A gas surface mass density is about Σ g ∼ 2 × 10 −4 − 2 × 10 −3 g cm −2 . The value of mass-to-flux ratio Λ is about 20 ∼ 2000 g cm −2 G −1 . Therefore, the λ parameter falls in the range of 0.03 ∼ 3. Here, we estimate the maximum range for every parameter in the galactic context while keeping in mind the constraint of 1 + ǫδ > 0 for the background magneto-rotational equilibrium. The variation range of ǫ may be as large as −1110 ∼ 0.9. Taking δ as 0.01, we choose the ǫ parameter to fall within the rough range of −100 ∼ 0.9. In a typical late-type spiral galaxy, the stellar velocity dispersion as is usually several times higher than the sound speed ag of the gaseous disc, such that β > 1 in our model analysis. We have already shown the Θ range to be 1 < Θ < 2. As Θ and β always appear together, we use a single combined parameter χ ≡ β/Θ and further take χ > 1 in reference to β > 1. As 1 < Θ < 2, the effect of Θ variation on χ is somewhat limited. For example, one may take ag = 7 km s −1 and as = 30 km s −1 for a typical late-type spiral galaxy. It is then justifiable to take χ = a 2 s /(Θa 2 g ) > 1. The introduction of B parameter by definition (113) is especially useful to estimate effects of an isopedic magnetic field. It is clear that inequality B ≥ 1 holds valid given the condition ǫ ≤ 1 and the requirement 1 + ǫδ > 0 for the background magneto-rotational equilibrium. With a fixed value of δ, the increase of the isopedic magnetic field flux will lead to an increase of B. It turns out that B is a fairly good indicator for the isopedic magnetic field, although the relationship between the two is not as simple as a linear one.
Formally, the two real y ≡ D 2 s solutions of quadratic equation (108) are given by We now introduce an important parameter Ra (see Appendix A for more details) defined by Based on the mathematical conclusion reached in Appendix A, for the y1 solution with the plus sign in (115), we have y1 > Ra > 0 when C2 > 0 and y1 < −1 when C2 < 0.
For the y2 solution with the minus sign in (115), we have −1 < y2 < Ra.
In the previous analyses , 2004aLou & Zou 2004), the two solutions y1 and y2 are sometimes referred to as supersonic and subsonic rotation solutions respectively, because as rotational Mach numbers, y1 is almost always greater than 1 and y2 remains always less than 1. In comparison, the situation becomes somewhat different in the presence of an isopedic magnetic field. In particular, y1 may drop below 1 for some specific values of m, δ, B and χ; this could also happen when there is no magnetic field. Furthermore, with the boosting of the magnetic field, y2 may become greater than 1 even under the condition χ ≥ 1. Note that y2 remains always less than 1 in the absence of magnetic field. Thus for proper terminologies, we use the F-wave solution for y1 and the S-wave solution for y2, because when they are both physical, y1 and y2 remain always faster (F) and slower (S), respectively.
When B = 1 for the absence of magnetic field, we have Θ = 1 and hence χ = β. This situation is exactly the same as the problem analyzed by  which will be frequently referred to in our following discussions.
For the F-wave solution y1, there exists an upper limit for the magnetic field intensity. When the magnetic flux becomes so strong such that B has a sufficiently high value, C2 will become negative and thus y1 becomes unphysical for being negative. The physical constraint on B for the F-wave solution y1 is therefore where the lower bound on the left-hand side (LHS) is from the definition of B while the upper bound on the right-hand side (RHS) is derived by requiring C2 to be positive. We note that with the increase of both m and δ, the allowed range of B for the existence of an F-wave solution y1 becomes enlarged. Moreover for y1 > 0, y1 decreases with increasing χ.
For the S-wave solution y2, the somewhat loose constraint is that −1 < y2 < Ra. For a physical solution, we should require y2 > 0. Since y1 > 0 for C2 > 0 and y1 < 0 for C2 < 0, we infer from the solution property y1y2 = C0/C2 that y2 > 0 for C0 > 0 and y2 < 0 for C0 < 0, respectively. For a physical S-wave solution y2 > 0, the condition C0 > 0 is simply equivalent to the following inequality When B = 1 for the absence of magnetic field, by setting the LHS of the above inequality (118) equal to zero, we determine a critical value χc for χ which is the same as βc parameter of ; as β or χ is increased to become greater than βc or χc, the S-wave solution y2 will change from positive to negative values. Furthermore, the increase of B leads to an increase of χc. For B > 1, the critical value χc is given explicitly by One interesting point to emphasize is that when B becomes sufficiently large, there is no limit on χ and the S-solution y2 remains always positive. More specifically, when the following condition (120) is fulfilled inequality (118) will be always satisfied and thus the Ssolution y2 remains always positive. Without magnetic field, there always exists at least one positive solution (see . In contrast, when an isopedic magnetic field is anchored in a gaseous SID, there does exist a combination of parameters such that both y1 and y2 solutions may become negative. Under these circumstances, there would be no stationary MHD density wave patterns supported by SID rotation. For this to happen, the range of B parameter is given below (117) and inequality (120)]; the LHS of inequality (121) leads to y1 < 0, while the RHS of inequality (121) makes y2 < 0 possible and the χ parameter  We take δ = 0.2 (solid lines), 1 (dotted lines) and 5 (dash-dotted lines), respectively. The allowed range for χ is [1, +∞), and we plot these curves within the χ interval (0, 20]. This figure is the same as figure 1 of  and serves as our reference to understand the role and effects of an isopedic magnetic field.
For the phase relationship between the two surface mass perturbations µ s and µ g , we note that equation (102) can be written in the form of .
According to expression (122) of µ g /µ s , it is fairly clear that for the F-wave solution (i.e., y1 > Ra with the unphysical case of y1 < 0 being ignored), the ratio µ g /µ s remains always negative, while for the S-wave solution y2, the ratio µ g /µ s remains always positive. Physically, when the two density perturbations are out of phase, the gravity effect is weaker and the MHD density wave speed is faster, while when the two density perturbations are in phase, the gravity effect is stronger and the MHD density wave speed is slower (Lou & Fan 1998b).
In addition to the requirement of D 2 s > 0, we must also make sure of D 2 g > 0 such that the relevant MHD density wave mode is physically plausible. We now write equation (107) in a simple form of Without magnetic field with B = 1 and under the assumption of χ = β > 1, it is easy to show that D 2 g remains always  s versus χ variation for the aligned case with |m| = 2 and B = 4 for the presence of an isopedic magnetic field. We take δ = 0.2 (solid lines), 1 (dotted lines) and 5 (dash-dotted lines), respectively. Again, we use B in the figure to correspond to B in the main text. The vertical line marks the constraint of χ > 1. Note that the F-wave solution y 1 with δ = 1 is not shown as it drops far below zero. In order to give an intuitive comparison between Fig. 1 and Fig. 2 here, we use the same scale for the reference frame. It should be noted that only a range of χ is permitted when the S-wave solution y 2 is positive due to the constraint of χs.
positive according to equation (123) Lou & Zou 2004). However, this may not be true in the presence of an isopedic magnetic field with B > 1 and D 2 g may become negative even for χ > 1. In other words, one more constraint of D 2 g > 0 should be checked carefully. On the basis of Appendix A, we here define a new critical value χs and only when χ > χs can inequality D 2 g > 0 be guaranteed. In the absence of magnetic field (i.e. B = 1), we readily see that χs < 1. So a more general condition for a physical solution would be χ > max{1, χs}. This limit of χs emerges as a constraint in the presence of an isopedic magnetic field (B > 1). With the increase of the magnetic flux (i.e. B grows), the value of χs increases accordingly. By comparing the expressions of χc and χs, it should be noted that wherever there exists an upper limit χc, this χc remains always larger than χs.
We now review several properties of the two solutions of y ≡ D 2 s as a summary. In the absence of magnetic field with B = 1 (see Figure 1), our analysis here is the same as those of . The F-wave solution y1 remains always positive and the S-wave solution y2 changes from positive to negative as χ becomes larger than χc. In the presence of an isopedic magnetic field with B > 1, we define a χs instead of χc. For a physical configuration solution, we then require χ > max{1, χs}. With the increase of the magnetic flux (i.e., an increase of B), the two solutions y1 and y2 become larger and the value of χc increases accordingly. However, as B exceeds a certain given value [see expression (117)], the F-wave solution y1 becomes negative and thus unphysical (see Fig. 2). In addition to the requirement 1 + ǫδ > 0 on the magnetic field for the background magneto-rotational equilibrium, a further condition (117) is needed for the existence of an F-wave solution y1. As χs becomes larger than 1, the requirement of χ > max{1, χs} means that only if χs < χ < χc, can the S-wave solution y2 be physically valid. Note that these constraints on χ do not apply for the F-wave solution y1. As the magnetic field strength is increased further [see equation (120)], the upper limit χc for χ will disappear and only the lower limit χs exists (see Fig. 3). It is generally true that when the two solutions y1 and y2 are both positive and physical, they both increase with the increase of B and decrease with the increase of χ.
Physically, for the F-wave solution with µ g /µ s < 0, the surface mass density perturbations of stellar and gaseous SIDs are completely out of phase to reduce the effect of gravity. Meanwhile, a faster azimuthal density wave speed corresponds to a faster SID rotation (larger D 2 s ) in order to maintain a stationary MHD configuration. For the S-wave solution with µ g /µ s > 0, the surface mass perturbations in gaseous and stellar discs are in phase with each other. As the effect of self-gravity is enhanced, the azimuthal density wave speed becomes slower corresponding to a slower SID rotation (smaller D 2 s ) in order to sustain a stationary MHD perturbation configuration.
The isopedic magnetic field plays an interesting role. When the magnetic field grows stronger as B becomes larger, the F-wave solution y1 and the S-wave solution y2 both become larger. According to expression (122) for the ratio µ g /µ s , it is easy to demonstrate that µ g /µ s will decrease with the increase of D 2 s for both y1 and y2 solutions. For the F-wave solution, the decrease of µ g /µ s means an increase of the absolute value of µ g /µ s , while for the S-wave solution, this ratio magnitude becomes smaller. A physical interpretation is that for a given composite SID system, an isopedic magnetic field will amplify the contrast of surface mass density perturbations in the two SIDs for an F-wave mode and will reduce such contrast of surface mass density perturbations in the two SIDs for a S-wave mode.
As noted earlier Lou & Fan 2002;Lou & Zou 2004), we here emphasize again the physical perspective that stationary aligned MHD perturbation configurations should be regarded as purely azimuthal propagation of MHD density waves counterbalanced by the advection of (M)SID differential rotation. Together with equation (107), we may write equation (105) in the following form to derive the two y solutions.
where C2 and Ra are defined by expressions (110) and (116), respectively. For the F-wave solution y1, we have D 2 s > Ra when D 2 s > 0 which is guaranteed by C2 > 0. Moreover from equations (125) and (123), we can show that for a positive F-wave solution y1, the corresponding D 2 g is also positive such that the overall solution is physically plausible.
We now write solution condition (105) for stationary MHD configurations of a composite SID system in a physically more suggestive form of The RHS of equation (126) represents the mutual gravitational coupling between the 'fluid' stellar SID and the 'magnetofluid' of a gaseous isopedic MSID. Without this gravitational coupling, the two factors on the LHS would emerge as two separate solution conditions for stationary aligned perturbation configurations with |m| ≥ 2 of stellar SID and of gaseous isopedic MSID, respectively. For the stellar SID alone, the condition would be and for the gaseous MSID alone with an isopedic magnetic field, the condition would be m 2 Ω 2 g = κ 2 g + m 2 Θa 2 g /r 2 − 2πǫGΣ g 0 |m|/r .
According to the well-known dispersion relation of density waves first derived under the tight-winding or WKBJ approximation [Lin & Shu 1964, 1966Lin 1987; or equation (39) of Shu et al. 2000], dispersion relations (127) and (128) can be readily recovered by replacing the radial wavenumber |k| with the azimuthal wavenumber |m|/r and setting ω = 0 in an inertial frame of reference as noted by  in the study of stationary MHD perturbation configurations of a single MSID with a coplanar magnetic field. Along this line of reasoning, we know that dispersion relations (127) and (128) separately describe azimuthal propagations of hydrodynamic density waves in a stellar SID and of MHD density waves in a gaseous MSID with an isopedic magnetic field.

Secular Barlike Instabilities
We have constructed analytically stationary configurations for aligned nonaxisymmetric MHD perturbations in a com-  for the onsets of bar-type instabilities. Judging similarities and differences between a single (M)SID and a composite (M)SID system, we examine stability properties for stationary configurations of aligned and unaligned perturbations with an isopedic magnetic field in reference to the analysis of Shu et al. (2000).
The criteria for (secular and dynamic) barlike instabilities were sometimes expressed in terms of the ratio of the rotational kinetic energy T to the absolute value of the gravitational potential energy W (Ostriker & Peebles 1973;Binney & Tremaine 1987;Shu et al. 2000;Lou & Zou 2004). To obtain the scalar virial theorem in the present case of a composite MSID system, we start from the radial force-balance equations of the equilibrium state, where the axisymmetric dark matter halo potential is not included yet. An addition of equations (129) and (130) gives Multiplying equation (131) by a ring area element 2πr 2 dr and integrating from 0 to a finite radius L, we obtain where Here T is the rotational kinetic energy in the composite SID system, U is the equivalent 'thermal energy' contained in the composite MSID system including the effect of magnetic pressure, and W is the equivalent gravitational-work integral including the effect of magnetic tension. Using δ ≡ Σ g 0 /Σ s 0 , Ωs = asDs/r, Ωg = Θ 1 2 agDg/r together with equations (27) and (28), we cast these integrals in the following forms of For a composite MSID system of an infinite radial extent, all three integrals (136)−(138) above diverge as L → +∞ but their mutual ratios remain finite. For example, the ratio of the kinetic energy of disc rotation to the absolute value of the gravitational potential energy is This ratio can also be arranged into the following form which is explicitly symmetrized with respect to physical parameters of the two SIDs. Here, ǫ does not appear explicitly in the final expression. Note that the ratio T /|W| falls between 0 and 0.5 as usual (Binney & Tremaine 1987;Lou & Zou 2004) and increases with increasing D 2 s . For stationary configurations of aligned MHD perturbations in a composite system of two coupled SIDs, the two possible values of D 2 s give rise to two different values of T /|W| ratio; the larger and smaller values of D 2 s correspond to larger and smaller values of the T /|W| ratio, respectively. Furthermore, as the magnetic flux increases (i.e., an increase of B), the T /|W| ratio decreases. This simply means that the effect of an isopedic magnetic field tends to make the composite SID system more stable , 2004aLou & Zou 2004.
On the basis of an extensive numerical exploration for m = 2 to m = 10 and so forth, one can show the existence of a dividing line in the ratio T /|W| for global sta-  tionary configurations of aligned MHD perturbations. More specifically, for F-wave solutions y1, the ratio T /|W| ranges from m/(4m + 4) to 1/2, while for S-wave solutions y2, the ratio T /|W| ranges from 0 to m/(4m + 4) when m ≥ 2. Based on N −body numerical simulation experiments involving only ∼ 300 particles for a stability analysis of a rotating disc under self-gravity, Ostriker & Peebles (1973) suggested an empirical criterion that T /|W| ≤ 0.14 ± 0.02 is necessary but not sufficient against bar-type instabilities (see also earlier numerical simulation results of Miller et al. 1970;Hohl 1971). When T /|W| ≥ 0.14 ± 0.02, a disc system would evolve rapidly into bar-type configurations (e.g. Binney & Tremaine 1987). In the current problem of a composite MSID system, the lowest value of T /|W| ratio for F-wave solutions y1 is ∼ 0.1667, while for S-wave solutions y2, the T /|W| ratio can be lower than 0.14 ± 0.02.
In the above analysis, we have not included the effect of a massive dark matter halo. One obvious consequence of a massive dark matter halo is to increase |W| and thus decrease T /|W| ratio, and for both classes of global perturbation configurations, this implies the tendency towards a stability. At any rate, a comparison of the above two criteria suggests that given a background dark matter halo, F-wave configurations tend to be more stable than S-wave configurations.
By numerical experiments for m = 2 up to m = 10 and so forth, we found that when B = χ, the value of T /|W| ratio for S-wave solutions y2 is always m/(4m + 4).

Unaligned Logarithmic Spiral Disturbances
For global stationary configurations of unaligned or logarithmic spiral MHD disturbances, we take on the following set of exact density-potential pair relations that satisfy the Poisson integrals, namely µ l = σ l r −3/2 exp(iα ln r) , where σ l and v l are small constants for perturbations, α is a parameter to characterize the radial variation; the two small coefficients v l and σ l are algebraically related by with Nm(α) ≡ K(α, m) being the Kalnajs function (Kalnajs 1971) and superscript l = s, g. For logarithmic spiral perturbations in a composite MSID system, the above description is a sensible extension of earlier model studies. The radial scaling parameter α (closely related to the radial wavenumber) is naturally taken to be the same for perturbations in both stellar SID and gaseous isopedic MSID. In the following analysis and computations, we use two useful formulae of Nm(α). One is the recursion relation in m of Nm(α) for a fixed α value (Kalnajs 1971) and the other is the asymptotic expression of Nm(α) Nm(α) ≈ (m 2 + α 2 + 1/4) −1/2 in the regime of m 2 + α 2 ≫ 1 (e.g., Shu et al. 2000). When the requirement for accuracy is not so stringent in some quantitative analyses, expression (145) may even be used to compute logarithmic spiral solutions with |m| = 1. For |m| > 0, we now proceed to solve equations (59) and (60) that can be cast into more compact forms of where the four coefficients H1, H2, G1 and G2 are Substituting µ s and µ g in the form of expression (141) into equations (146) and (147), we readily obtain Multiplying both sides of the above two equations and re-moving µ s µ g = 0, we derive the solution condition as [1 − H1(m 2 + α 2 + 1/4)][1 − H2(m 2 + α 2 + 1/4)] = G1G2(m 2 + α 2 + 1/4) 2 . (154) Using definitions (148)−(151), condition (154) then becomes This equation may be rearranged into the form of Using expression (107) for D 2 g , we readily obtain a quadratic equation in terms of y ≡ D 2 s , where the three coefficients C2, C1 and C0 are In the following, we shall remove some nonzero common factors of C2, C1 and C0 to simplify these expressions. For this purpose, we introduce a new parameter Together with parameter B defined by equation (113) as in the aligned case, we obtain Consequently, we have the determinant ∆ of quadratic equation (155) above as Formally, there are two different real solutions y1 and y2 of quadratic equation (155) in the forms of for ∆ > 0. Based on our earlier discussion for the axisymmetric situation of m = 0 in equations (88) and (89) for both aligned and unaligned cases, the above results remain valid for the axisymmetric case of m = 0 in the unaligned case (i.e. perturbations with radial oscillations). Note that ∆ remains always positive for m ≥ 1. The rare case of ∆ = 0 can only happen when m = 0 and M + 1 = 0. In general, there are therefore two different real solutions y1 and y2 and these two solutions may become the same only if m = 0 and M+1 = 0.
We now proceed to discuss the three situations m ≥ 2, m = 1 and m = 0 in three separate subsections below.
As in the aligned case, we introduce a similar Rs parameter (see Appendix B for more details). On the basis of the mathematical analysis in Appendix B, for the y1 solution, we have y1 > Rs > 0 when C2 > 0 while y1 < −1 when C2 < 0, while for the y2 solution, we have −1 < y2 < Rs independent of the sign of C2. Here, the role of this Rs parameter for the unaligned case is similar to the role of Ra parameter for the aligned case. Note that when m ≥ 2, parameter Rs remains always positive. As in the aligned case, we still refer to y1 as the F-wave solution and y2 as the S-wave solution for the unaligned case, respectively.
To analyze the phase relationship between the two surface mass perturbations µ s and µ g in the coupled stellar SID and isopedic MSID, we recast equation (152) in the form of From expression (160) of ratio µ g /µ s , it is clear that for the F-wave solution y1 (i.e., y1 > Rs with the unphysical situation of y1 < 0 being ignored), ratio µ g /µ s is always negative, while for the S-wave solution y2, ratio µ g /µ s is always positive. Note that these results are qualitatively the same as those in the aligned case with m ≥ 2. The same as the aligned case, for the F-wave solution y1, ratio µ g /µ s < 0 means that surface mass density perturbations in stellar SID and gaseous MSID are out of phase while for the S-wave solution y2, ratio µ g /µ s > 0 means that the surface mass density disturbance in the gaseous MSID is in phase with the surface mass density disturbance in the stellar SID. It becomes transparent that in-phase and out-of-phase surface mass density perturbations are intrinsic characters for the two solutions y1 and y2 when m ≥ 2. Physically, both aligned and unaligned perturbations represent MHD density wave propagations in a composite MSID system. As MHD density wave patterns with out-of-phase surface mass density perturbations travel faster, a higher rotational Mach number Ds is required in order to strike a global stationary configuration. In contrast, MHD density wave patterns with in-phase surface mass density perturbations travel slower, a lower rotational Mach number Ds is needed in order to strike a global stationary configuration. This phase relationship for the two surface mass density perturbations is retained even in the presence of an isopedic magnetic field (Lou & Fan 1998;Lou & Zou 2004). The solution y stands for D 2 s which is the square of the rotational Mach number Ds of the stellar SID. The two solutions y1 and y2 correspond to the square of two different rotational Mach number. One obvious physical requirement is that both y1 and y2 must be positive. Meanwhile, we also need to require D 2 g to be positive in order to establish a plausible composite MSID system. According to equation (123), we have the correspondence between the condition D 2 g > 0 and the condition We discuss in the following these two constraints for F-wave solutions y1 and S-wave solutions y2, separately. The two constraints must be satisfied simultaneously to limit the parameter space, otherwise there are no physical y solutions. For F-wave solutions y1 to exist, we have the following mathematical conclusions from Appendix B, namely are guaranteed by C2 > 0. As a result, we can set a limit on the strength of an isopedic magnetic field by requiring C2 > 0. For the B parameter, we thus have inequalities for the existence of a physical F-wave solution y1. As in the aligned case, this constraint means that too strong an isopedic magnetic field would be impossible for sustaining a stationary F-wave solution y1. By increasing both m and δ and decreasing α, the allowed range of B for sustaining F-wave solution y1 becomes enlarged. Once the F-wave solution y1 exists, we always have a decreasing y1 as both χ and B increase (see Appendix B for details). Through our analysis, we also find that for B > 1 and α → ∞, the C2 coefficient will become negative and thus lead to a negative y1. Therefore, there exists a critical value αc of α such that y1 > 0 exists only when α < αc; this critical αc is For the S-wave solution y2 > 0 to exist, we establish the following two correspondences between y2 > 0 and C0 > 0 , and between y2 > B χ − 1 and respectively. The requirement of C0 > 0 for a positive y2 leads to a constraint on χ parameter, that is, χ < χc where the critical value χc for χ is explicitly defined by We further note that when B parameter becomes sufficiently large to satisfy the following inequality the χ parameter will be free from the constraint of χ < χc because such χc no longer exists. From the inequality we derive another constraint on χ, that is, χ > χs where χs is explicitly defined by We can readily demonstrate that χs < χc such that physical s versus α of the unaligned case with specified δ = 0.2, χ = 1.5 and B = 1.5, the two solid solution curves y 1 and y 2 are for m = 2 and the one dash-dotted solution curve y 2 is for m = 1. Here, we use B in the figure to correspond to B in the main text. All these solutions increase with increasing α. The dashed vertical line represents α = αc. When α becomes greater than αc, the branch of y 1 (m = 2) solution flips and becomes negative. In comparison, we see that y 2 with m = 1 is greater than y 2 with m = 2. solution y2 > 0 could exist for χs < χ < χc. However, we show presently that χs > χc may happen in certain situations and the two requirements are incompatible with each other; there is thus no physical solution for a positive y2.
The allowed range for χ parameter is χ > 1. When B = 1 for the absence of an isopedic magnetic field, the critical χs remains always less than 1 such that the second requirement will be satisfied automatically. As the strength of an isopedic magnetic field gets stronger, the χ > χs constraint becomes indispensable for a positive y2.
Our analysis above bears a strong resemblance to that for the aligned case. In fact, should we set the parameter α = 0 for a purely azimuthal propagation of MHD density waves, we can approximately obtain nearly the same results as those of the aligned case. This naturally suggests that the aligned case should be only regarded as a special case of the unaligned cases Lou & Fan 2002;Lou & Zou 2004).
Through extensive numerical explorations, we can show empirically that both the F-wave solution y1 and the S-wave solution y2 increase with increasing α (see Figure 5). From the perspective of global stationary MHD density waves, this is a physically sensible result.
By numerical experiments (see Appendix B), we have C2 < 0 and only y2 solution being positive. The ratio µ g /µ s of y2 solution remains always positive with in-phase surface mass density perturbations. We also know that y2 increases with decreasing χ and increasing B. Numerically, we find that y2 solution increases with increasing α. The two constraints for sensible y2 solutions may be expressed as respectively (see Appendix B2). The mathematical forms of these constraints are too involved to be shown here. Note that properties of y2 solution with m = 1 are just like those of y2 in the cases of m ≥ 2 (see Fig. 6). As the y2 solutions of m = 2 and m = 1 are fairly similar to each other, the y2 solution changes smoothly as m increases sequentially.
On the basis of extensive numerical experiments, we found that the five expressions below will become negative when α is sufficiently small. As α is increased across several different critical values, they will in turn become positive. We now proceed to discuss consequences of this α variation in different ranges separately. The analysis is fairly complicated and we summarize the basic results here. By imposing the two constraints of y > 0 and y > (B/χ) − 1, we found that when α < 1.113 and α > 1.793, there is only one y solution that is physically valid under certain conditions. For 1.113 ≤ α ≤ 1.793, there is no physically valid solution at all. We then find that the valid y solution is characterized by in-phase surface mass density perturbations (i.e. a positive µ g /µ s ratio). As before, the y solution increases with decreasing χ and increasing B in both ranges of α < 1.113 and α > 1.793. Numerically, we again find that when α < 1.113, the solution of D 2 s decreases with increasing α. In comparison, when α > 1.793, D 2 s first decreases with increasing α and upon crossing a certain point, D 2 s increases with increasing α. For the positive portions of solution D 2 s without going into the extreme, the basic profile of D 2 s versus α is qualitatively similar to the results of the single SID in Shu et al. (2000) and of the composite SIDs in . In fact, our model analysis is more general to include the results of Shu et al. (2000) and of .
From the dispersion relation of MHD density waves in the WKBJ or tight-winding approximation , 2004bLou & Zou 2004), we see that the solution curve in the D 2 s versus α profile represents the marginal stability curve that separates the regime of stable oscillations from the regime of unstable oscillations. As noted by Shu et al. (2000), when α < 1.113, the parameter regime under the curve in the lower-left corner represents the rotation and magnetic field modified Jeans collapse regime. As α stands for the radial wavenumber and m stands for the azimuthal wavenumber, a small α corresponds to a long-radial wavelength MHD density wave disturbance. When α > 1.793, the parameter regime above the solution curve represents the ring fragmentation regime for large α corresponding to short-radial wavelength MHD density wave (see Fig. 7).
As χ parameter increases, the y solution decreases. Hence the collapse regime shrinks and the ring fragmentation regime grows with increasing χ (see Fig. 8). As B parameter increases, the y solution increases. Hence the collapse regime grows and the ring fragmentation regime shrinks with increasing B (see Fig. 9).
Physically, χ parameter stands for the ratio in the effective temperatures of the two SIDs and B represents the strength of the isopedic magnetic field. The above results mean that when the 'fluid' stellar SID becomes much 'hotter' than the gaseous isopedic MSID, the stellar SID tends to be more stable against collapse for large-scale disturbances but less stable against ring-like fragmentation for small-scale disturbances. For a stronger isopedic magnetic field, the tendency of collapse becomes easier while the ring-fragmentation becomes more difficult. Here the B parameter involves the background equilibrium parameters. While we seem to vary only the background magnetic field strength, the background magneto-rotational equilibrium actually changes accordingly.

CONCLUSIONS, SUMMARY AND DISCUSSIONS
For modelling an actual spiral galaxy, our treatment of a composite system of MSIDs, while highly idealized still, contains several physically realistic key ingredients. The formulation and results of our analysis provide the basic rationale for developing useful astrophysical concepts. The two theorems of Shu & Li (1997) have now been generalized to a composite MSID system with an isopedic magnetic field. Instead of an assumption, the isopedic relation is shown here to be an initial condition maintained by ideal MHD equations. In particular, we have taken into account of long-range effects of gravity in full without the usual local WKBJ approximation. We have succeeded in constructing various global perturbation configurations, such as bars, spirals, barred spirals etc. By analogies, we proposed specific transition criteria among different stationary perturbation configurations.
With the isopedic magnetic field geometry, we have set the stage to further study global patterns of circumnuclear and spiral arm galactic winds. We now discuss plausible consequences and possible applications to spiral galaxies from several perspectives. On much smaller scales, physical processes involve various star formation activities and their interactions with the surrounding ISM. As luminous massive OB star formation tends to be more numerous in regions of relatively high gas concentration and of stronger magnetic field, the collection of numerous young bright stars will then outline a spiral pattern embedded in relatively high density spiral arms. Theoretical studies on the large-scale dynamics of galactic spiral structure are closely tied to a chain of star formation processes and to the physical manifestation of spiral MHD density waves. For the global star formation rate in a disc galaxy (e.g., Kennicutt 1989;Elmegreen 1994;Silk 1997), magnetic field plays an important role in the gaseous disc in the context of large-scale instabilities (e.g., Lou & Fan 2000;Lou et al. 2001). We would emphasize that should largescale galactic disc instabilities be indeed responsible for the global star formation rate, then analyses and simulations of a composite disc system involving magnetic field would be more realistic and physically meaningful. The oft-used axisymmetric stability criterion involving the Q parameter is significantly modified in a composite disc system with an isopedic magnetic field Lou & Zou 2004).
Galaxies form and evolve on large scales. To provide a dynamic basis for the classification of galactic morphologies is one major goal of galactic research (e.g., Bertin et al. 1989a, b;Bertin & Lin 1996). The modal approach to this problem has gained remarkable progress in explaining various configurations of spiral, barred-spiral galaxies with observational constraints. For example, a barred galaxy is expected to be associated with a relatively heavy disc whose halo mass inside the optical radius is small. The potential of a rotating bar can give rise to quasi periodic perturbations in the disc and hence affect the entire configuration of a disc. While limited by the stationarity requirement in our combined analytical and numerical solution procedure, we are able to take into account of long-range effects of gravitation in full for MHD perturbations in a composite MSID system without the deficiency of the usual local WKBJ approximation. For example for m = 2, we have constructed separately barred configurations (aligned) and logarithmic spiral configurations (unaligned). By the principle of linear superposition, we might be able to also construct barred logarithmic spiral configurations in a composite MSID system by choosing appropriate sets of parameters. In addition to theoretically modelling galactic configurations, these solu-tions are valuable in testing, initializing and benchmarking hydrodynamic and MHD numerical simulation codes for further nonlinear studies.
As important extension and generalization of the earlier theoretical analyses of Shu et al. (2000) on zero-frequency (i.e. stationary) aligned and unaligned MHD perturbation configurations of an isopedically magnetized SID and of  on stationary aligned and unaligned perturbation structures in a composite system of two-fluid SIDs, we have constructed analytically stationary configurations of aligned and unaligned MHD perturbations in a composite system of a 'fluid' stellar SID and a 'magnetofluid' gaseous MSID. Other closely relevant work include Galli et al. (2001), , Lou & Fan (2002), , Shen & Lou (2004a, b), Chakrabarti et al. (2003), Shen, Liu & Lou (2005), Lou & Zou (2004. While this composite model of one SID and one MSID is highly idealized in many aspects, it does contain several necessary and more realistic elements that are crucial to understand large-scale structures and dynamics of spiral galaxies. From the basic fluid-magnetofluid equations for a composite system of a stellar SID and a gaseous MSID, we begin with an axisymmetric background in a magneto-rotational equilibrium. The necessary condition of 1 + ǫδ > 0 arises here in order to maintain an equilibrium with a positive gas surface mass density Σ g (eqn. 28) which in turn puts an upper limit on the magnetic flux under the isopedic condition of constant Σ g /Bz. One important conclusion of our analysis is that a constant Σ g /Bz is an initial condition sustained by ideal MHD equations during disc evolution. Physically, the magnetic tension force cannot be too strong to make an equilibrium of the gaseous MSID impossible (eqn 26) in the first place. Massive dark-matter halos tend to stabilize disclike galaxies as expected (Ostriker & Peeble 1973;Hohl 1976;Miller 1978). Although quantitative results for a partial composite MSID system are not detailed here, we suggest that the potential of a dark matter halo will certainly stabilize a composite MSID system and allow for a stronger flux of isopedic magnetic field based on our working experience. The B parameter [see definition (113)] contains this constraint explicitly.
For MHD perturbations, we derive two sets of linearized equations in the stellar SID and the gaseous MSID. Setting ω = 0, we derive the stationary MHD dispersion relation with the gravitational coupling. By properly choosing different potential-density pairs to satisfy the Poisson integral, we examine the two classes of aligned and unaligned logarithmic spiral MHD configurations in parallel. As a result of the gravitational coupling between the two SIDs, we obtain a quadratic equation in terms of y ≡ D 2 s to construct stationary MHD perturbation configurations of the composite MSID system, namely, equation (108) for the aligned case and equation (155) for the unaligned logarithmic spiral case. The two necessary physical requirements are that both D 2 s and D 2 g should be positive. On the basis of rigorous mathematical derivation and analysis, we have reached several conclusions and results.
Meanwhile, we have explored parameter regimes numerically to reveal a few trends of parameter variations empirically.
Let us first summarize the basic results of the aligned case. For axisymmetric MHD perturbations with m = 0, the resulting configuration is simply a rescaling of one axisymmetric equilibrium to a neighbouring axisymmetric equilibrium which is certainly allowed but trivial for our purpose. For nonaxisymmetric MHD perturbations with m = 1, the resulting quadratic equation of y ≡ D 2 s is automatically satisfied for arbitrary values of y (Shu et al. 2000;. For nonaxisymmetric MHD perturbations with m ≥ 2, there are two sets of y ≡ D 2 s solutions representing two different classes of stationary global MHD perturbation configurations. For the phase relationships between the surface mass density perturbations in the two SIDs, the lower Ds solutions are featured by an in-phase relation while the higher Ds solutions are featured by an out-of-phase relation. For a weak isopedic magnetic field, the out-of-phase solution always exists (D 2 s > 0 and D 2 g > 0 simultaneously), while for the existence of the in-phase solution, the condition χs < χ < χc must be met with χs and χc defined by equations (124) and (119) for the aligned case and by equations (163) and (165) for the unaligned case. As χs < 1 may happen as the isopedic magnetic field becomes weak, the requirement χ > χs will be automatically satisfied by the physical constraint χ > 1. As the isopedic magnetic flux increases across a certain value, the out-of-phase solution vanishes while the constraint of the in-phase solution becomes χ > χs with χs > 1 and with the constraint χ < χc being irrelevant. Moreover, we emphasize that in order to have D 2 g > 0, the condition χs < χ must be satisfied for both aligned and unaligned cases and for both F-wave and S-wave solutions; the constraint χs < χ for F-wave solutions is satisfied automatically once other requirements are met.
In the framework of the density wave theory, the case of aligned perturbations represents a purely azimuthal propagation of density wave  and the stationarity is sustained by the advection of a counter SID rotation. In general, the class of unaligned logarithmic spiral perturbations represents MHD density waves with both radial and azimuthal propagations where parameter α is an effective radial wavenumber. Another special situation is a purely radial propagation with m = 0 in the 'unaligned' class. It is thus not surprising that with α = 0 in the unaligned class, we obtain fairly similar results in parallel with the aligned class for m ≥ 2. For example, by setting α = 0 in the unaligned class, we arrive at essentially the same results as in the aligned class for m ≥ 2. This justifies our physical interpretation that the aligned class is merely a special case to the unaligned class. For the unaligned class, we obtain two solutions of D 2 s characterized by out-of-phase and in-phase surface mass density perturbations, respectively. Other solution properties of the unaligned class remain qualitatively the same as those of the aligned class.
The unaligned cases with m = 1 and m = 0 need more special considerations. For m = 1, only the D 2 s solution with in-phase surface mass density perturbations could be valid under certain constraints. For m = 0, only the D 2 s solution with in-phase surface mass density perturbations could still exist; when 1.113 ≤ α ≤ 1.793, there is no physical solution though. From the local WKBJ dispersion relation, the profile D 2 s versus α represents the marginal stability curve , 2004aLou & Zou 2005;Shen, Liu & Lou 2005). In general, when D 2 s and α are both sufficiently small, rotation and magnetic field modified gravitational instabilities for Jeans collapse occur and when D 2 s and α are both sufficiently large, MHD ring fragmentation instabilities appear. Between these two unstable regimes is the stable regime of SID rotation for axisymmetric oscillations. By analytical analyses, we find that all physically valid D 2 s solutions increase with decreasing χ and increasing isopedic magnetic flux. This is a quite general conclusion for both aligned and unaligned cases. The substantial mathematical procedures for these conclusions can be found in Appendices A and B.
Finally, we discuss several possible implications and applications of our model analysis. After the theoretical prediction of the solar wind by Parker (1958Parker ( , 1963 and the subsequent confirmation by spacecraft observations in the early 1960s, the more general concept of stellar winds and the relevant theories have been developed in parallel (e.g., Burke 1968;Holzer & Axford 1970;Lamers & Cassinelli 1999). Meanwhile, the concept of galactic winds has gradually emerged (e.g., Johnson & Axford 1971). Intrigued by the almost absence of ISM materials in most elliptical galaxies, Mathews & Baker (1971) also came to the initial concept of galactic winds. It is now widely believed that a galactic wind (sometimes also referred to as a galactic superwind) is powered by starburst activities and supernova explosions (e.g., Chevalier & Clegg 1985;Heckman et al. 1990). Because of various observational limitations, the detection of a galactic wind has remained a challenge for many years. Along with the advance of observational technologies, evidence for galactic winds have gradually surfaced in the multi-wavelength studies. Using the ROSAT data for Xray emissions, Strickland et al. (1997) presented a galactic wind picture for the nearby galaxy M82. Matthews & de Grijs (2004) found the evidence for a galactic wind in the edge-on Sbc galaxy UGC 10043 by optical imaging and spectroscopy. Melo et al. (2003) reported the detection of supergalactic winds in the edge-on starburst galaxy NGC 4631. Together with galactic-scale outflows, a relevant concept and phenomenon of superbubbles have been proposed and observed (e.g., Tomisaka, Ikeuchi & Habe 1981;Ferrière 2001). Both phenomena are the results of supernovae, hypernovae, violent burst activities of massive star formation in the densest regions of host galaxies. Galactic winds are important processes to channel ISM materials along open magnetic fields into the intergalactic medium, while superbubbles fail to get out of the reign of the host galaxies (e.g., Tenorio-Tagle, Silich & Muñoz-Tuñón 2003). The physical reasons why superbubbles cannot become galactic winds remain unclear.
Based on our model analysis, we here propose that the global structure and geometry of galactic magnetic field should be mainly responsible for the coexistence of largescale outflows and superbubbles in disc galaxies. In strong analogy to similar phenomena such as global coronal mass ejections, violent solar flares from active regions with closed magnetic field over the solar surface and fast solar wind streams coming out of solar coronal hole regions with open magnetic fields, we suspect that bubbles and superbubbles may involve closed magnetic fields anchored at the galactic disc plane, while large-scale galactic winds pumped by energetic sources must be guided by open magnetic field lines to reach outer regions. The model analysis here and those of Lou & Zou (2004 and of Shen, Liu & Lou (2005) are complementary to each other in terms of overall magnetic field geometries associated with discs. In the case of a coplanar magnetic field, Parker instabilities inflated by galactic cosmic rays and supernova explosions can lead to formation of bubbles and superbubbles. In the case of an isopedic magnetic field, Velikhov-Chandrasekhar-Balbus-Hawley instabilities, starburst activities and the cosmic-ray pressure together pump galactic outflows from circumnuclear regions and along spiral arms. We have separately prepared MHD stages in an idealized manner for further research on pertinent physical processes. In a real disc galaxy, these two different magnetic field geometries are intermingled and randomly distributed all over the disc. Intuitively, it is natural to imagine that outflows will be stronger at places where open magnetic field flux is high. As large-scale dense regions in a galactic plane, the spiral arms in a galaxy are expected to carry more magnetic field flux because of the frozen-in condition for magnetic flux in the gas and of small-scale MHD dynamo processes. Comparing other regions in the disc (except for the nucleus and circumnuclear regions), stronger outflows are expected to emerge from spiral arm regions. We refer to this scenario as Spiral Arm Galactic Winds.

APPENDIX A:
In the aligned cases of m ≥ 2, we defined two parameters In the quadratic equation of y ≡ D 2 s for a stationary global MHD perturbation configuration the three coefficients C2, C1 and C0 are defined by From quadratic equation (A3), one can readily derive the expression of the determinant ∆ and show that There are thus two real y ≡ D 2 s solutions for equation (A3). For the convenience of analysis, we introduce a useful and important Ra parameter, The two following inequalities hold, namely According to the Viete theorem (e.g., John & Horst 1998) for a quadratic equation, we know that y1 > Ra > y2 > −1 for C2 > 0, while Ra > y2 > −1 > y1 for C2 < 0. Note that when C2 < 0 (i.e., X < 0), one must also have C1 < 0. It is then clear that y1 remains positive and negative for C2 > 0 and C2 < 0, respectively. To determine the phase relationship between the two surface mass density perturbations µ g and µ s , we derive where H1 and G1 are separately defined by expressions (98) and (100). Here, we are primarily interested in a y solution which is at least positive; a negative y solution is unphysical. As the case of y1 > 0 always corresponds to y1 > Ra > 0, it follows that µ g /µ s corresponding to y1 solution remains always negative for µ g and µ s being out of phase. As y2 < Ra for any values of C2, it follows immediately that the ratio of µ g /µ s corresponding to y2 solution remains always positive for µ g and µ s being in phase.
In addition to the physical requirement of y ≡ D 2 s > 0, we must also impose D 2 g > 0. In the purely hydrodynamic case of  without an isopedic magnetic field as studied here, D 2 g > 0 follows from D 2 s > 0 automatically. According to the background magneto-rotational equilibrium condition of radial force balance, we readily derive the inequality In other words, besides D 2 s > 0, we should also require D 2 s > B/χ − 1 in the case of a positive RHS.
We now explore the relations between the two real y solutions and their variations with specific parameters such as χ and B. For a chosen parameter, we apply the procedure of taking the derivative of the relevant y solution with respect to this specific parameter, namely where the prime ′ indicates a derivative with respect to this chosen parameter. It follows immediately that where y stands for either y1 or y2 solution given by By simple manipulations, we immediately arrive at where ∆ is the determinant of quadratic equation (A3). By these expressions, we examine the relations between the two y solutions and their dependence on relevant parameters. We first discuss solution properties of y1 and y2 by varying parameter χ and thus have explicitly It follows immediately that By expressions (A8) and (A9), the first derivatives y ′ 1 and y ′ 2 then take the forms of On the basis of our analysis following equation (A4), we have y ′ 1 < 0 and y ′ 2 < 0 for C2 > 0, and y ′ 1 > 0 and y ′ 2 < 0 for C2 < 0. The basic conclusion is that when y1 and y2 are both positive and thus physical, we have both ∂y1 ∂χ < 0 and ∂y2 ∂χ < 0 .
We next examine the solution properties of y1 and y2 as parameter B varies and define explicitly It then follows immediately that .

By rearrangements and manipulations, it is clear that
where the three coefficients D2, D1 and D0 are defined by .
By the properties of the two real solutions y1 and y2 and the expression above, we can readily show that D2y 2 1 + D1y1 + D0 > 0 and D2y 2 2 + D1y2 + D0 < 0, respectively. With the expression below from the second line of equation (A12) we immediately conclude that We note that for a given value of δ, parameter B is a function of ǫ, while parameter χ ≡ β/Θ is yet another function of ǫ according to the definition of Θ. It is then possible to write Θ as a function of B in the explicit form of .
As we have Note that for C2 < 0 and y1 < 0, the case is not considered here because of an unphysical negative y = D 2 s . Finally, we proceed to derive analytically the relevant constraints of requiring D 2 s > 0 and D 2 g > 0, which are equivalent to y > 0 and y > (B/χ − 1), respectively.
We then consider y2. As we know y1 > B/χ − 1 (equivalent to D 2 g > 0) and the properties of a quadratic equation, we therefore have the following two correspondences separately. For the following two inequalities to be valid, we clearly need to require that and respectively.

APPENDIX B:
For the unaligned cases or the logarithmic spiral cases, we first introduce a new parameter M as already defined in the main text, namely M ≡ m 2 − 2 (m 2 + α 2 + 1/4) .
In the quadratic equation we have the following forms of C2, C1 and C0 coefficients The determinant ∆ of quadratic equation (B1) can then be explicitly shown as so that two real y solutions are guaranteed.
In parallel to what we have done for the aligned case in Appendix A, we would like to find out the relationships between the y solutions and their variations with relevant parameters. Again, we use y1 and y2 to represent the two real y solutions defined by Following the same procedure of Appendix A, we derive where the prime ′ denotes the first derivative of a quantity with respect to a specific parameter. We first discuss the y solution properties with respect to the variation of parameter χ and define explicitly It follows immediately that In parallel with the analysis of the aligned case, we here introduce a similar parameter Rs in the form of It is then straightforward to demonstrate that We now examine behaviours of y solutions by varying parameter B. Following the same procedure, we write and derive explicitly It is straightforward to demonstrate that The phase relationship between the two surface mass density perturbations µ g and µ s is then given by , where H1 and G1 are separately defined by (148) and (150). In order to have both D 2 s and D 2 g being positive, we still have the same two requirements of D 2 s = y > 0 and D 2 s = y > B/χ − 1 simultaneously as in the aligned case.
In the following, we analyze the three situations of m ≥ 2, m = 1, and m = 0, separately.
According to the basic theory of a quadratic equation, we infer that y1 > Rs > y2 > −1 for C2 > 0, while Rs > y2 > −1 > y1 for C2 < 0. This inference is entirely similar to that of the aligned case. And we can get the similar condition for y2 that y2 > 0 ⇔ C0 > 0. Note also that Rs increases with increasing δ and α. By these analyses, we reach the following conclusions. For y1 and y2 being positive and thus physical, we have ∂y1 The perturbation mass density ratio µ g /µ s corresponding to y1 remains negative with µ g and µ s being out of phase, while µ g /µ s corresponding to y2 remains positive with µ g and µ s being in phase. The result obtained and the procedure taken here are just the same as those for the aligned case. Moreover, by setting α = 0 in the unaligned case, the equation for the unaligned case can be simplified to the same equation of the aligned case, implying that the aligned case is just a special example of the unaligned case.
By straightforward algebraic manipulations, we find that M + N1(α) > 0 and M + 1 > 0. As B ≥ 1, C2 remains always less than 0 and there must be a negative y solution (in fact this y solution must also be less than −1) that corresponds to y1. In this case, only the y2 solution may be positive in order to be physical. Because of y1y2 = C0/C2 and the solution property of y1 and C2 derived earlier, the physical requirement of y2 > 0 is clearly equivalent to C0 > 0 .
One can readily show the equivalence between the two inequalities Rs > 0 and M + N1(α)/(1 + δ) > 0. When δ < 0.145307118, Rs remains always positive for arbitrary values of α. Furthermore, we know that Rs > y2 > −1 for Rs > 0, while y2 > −1 > Rs for Rs < 0. It then follows that ∂y2 ∂χ < 0 and ∂y2 ∂B > 0 ; and the ratio of µ g /µ s corresponding to y2 remains always positive with µ g and µ s being in phase. This phase relationship is just the same as that of y2 solution in the case of m ≥ 2. Here, we have again the following two equivalent inequalities between y2 > 0 and C0 > 0 , and between respectively.
We now proceed to discuss each α−range separately.
For the y1 solution, we only need to think of the requirement y1 > (B/χ) − 1. After some analysis, we can demonstrate the equivalence of two inequalities for (B/χ) − 1 > Rs. For the y2 solution, we found by analytical derivation that the two inequalities y2 > 0 and y2 > (B/χ) − 1 are incompatible with each other. Thus in this α−range, there is no physical y2 solution.
We note that Rs = 0 for the special case of α = α1, so that y1 > Rs = 0 > y2 > −1. Here, only the y1 solution is positive and thus physical with the ratio of µ g /µ s being positive for in-phase µ g and µ s .
In this case with C2 > 0 and −0.7444 < Rs < 0, we have y1 > Rs > y2 > −1, such that y2 is negative and unphysical. The surface mass density ratio of µ g /µ s of y1 solution is positive with µ g and µ s being in phase. We further infer ∂y1 ∂χ < 0 and ∂y1 ∂B > 0 .
Through an analytical analysis, we can demonstrate the equivalence of two pairs of inequalities between for (B/χ) − 1 > Rs, and between y1 > 0 (B11) and C0 < 0 .
For the special case of α = α2, we still have a positive y1 solution. However, the conditions y1 > 0 and y1 > (B/χ) − 1 are incompatible with each other and there is then no physical y solution.
In this case with C2 > 0 and −1 < Rs < −0.7444, we have y1 > Rs > y2 > −1. The y2 solution remains negative and thus unphysical. The perturbation surface mass density ratio of µ g /µ s for the y1 solution is positive with µ g and µ s being in phase. We further infer the following inequalities ∂y1 ∂χ < 0 and ∂y1 ∂B > 0 .
The two inequalities y1 > 0 and y1 > (B/χ) − 1 are incompatible and there is no physical y1 solution in this α−range.

B3.4 The Range of α3 < α < α4
In this α−range, we have the following inequalities With C2 > 0 and Rs < −1, we have −1 > y1 > Rs > y2. Both y1 and y2 solutions are negative and thus unphysical. For α = α4, we have C2 = 0, corresponding to the divergence point of the quadratic equation; and only one negative y solution exists.
The negative y1 solution is unphysical. The perturbation surface mass density ratio µ g /µ s of y2 solution is positive with µ g and µ s being in phase. We further infer ∂y2 ∂χ < 0 and ∂y2 ∂B > 0 .
By an analytical analysis, we can demonstrate the equivalence of the following two sets of inequalities between and between y2 > 0 (B15) and C0 > 0 , respectively. For α = α5, we have C2 < 0 and thus a negative y1 solution. Only the y2 solution has the possibility of being positive. The ratio µ g /µ s of y2 solution for α = α5 is positive with µ g and µ s being in phase.

B3.6 The Range of α > α5
In this α−range, we have the following five inequalities With C2 < 0 and Rs > 0, we have Rs > y2 > −1 > y1. The negative y1 solution is unphysical. The perturbation surface mass density ratio µ g /µ s of the y2 solution is positive with µ g and µ s being in phase. We further infer ∂y2 ∂χ < 0 and ∂y2 ∂B > 0 .
By an analytical analysis, we can demonstrate the equivalence of the following two sets of inequalities between and between y2 > 0 (B19) and C0 > 0 , respectively. Finally, we give a complete summary for the m = 0 case with different ranges of α values. Physically, we could only have one real positive y solution with surface mass density perturbations µ g and µ s being in phase. For 0 < α < α2 = 1.113 and some constraints, y1 is the only physical solution. For 1.113 = α2 ≤ α ≤ α4 = 1.793, no physical solution exists. For α > α4 = 1.793 and some constraints, y2 is the only physical solution. All the physical y solutions increase with increasing B and decreasing χ.