Instabilities in two-fluid magnetized media with inter-component drift

Abstract

We analyse the stability of a magnetized medium consisting of a neutral fluid and a fluid of charged particles, coupled to each other through a drag force and exposed to differential body forces (for example, as the result of radiation forces on one phase). We consider a uniform equilibrium and simple model input physics, but do not arbitrarily restrict the relative orientations of the magnetic field, slip velocity and wavevector of the disturbance. We find several instabilities and classify these in terms of wave resonances. We briefly apply our results to the structure of SiO maser regions appearing in the winds from late-type stars.

instabilities, MHD, stars: winds, outflows, ISM: clouds, dust, extinction, ISM: kinematics and dynamics

1 Introduction

Multiphase flows are a widespread and important phenomenon in astrophysics. The difference between heating and cooling rates in different components of astrophysical gases often leads to the formation of a multicomponent medium, in which several phases with widely separate temperatures coexist near to pressure equilibrium. Effective multiphase behaviour can also result from differential coupling of distinct particle species to local magnetic fields or radiation driving forces. To first order, these differential forces will lead to drift velocities between the different components of the fluid, limited by the effect of frictional terms. However, this means that there is local source of free energy in the flow.

Radiation pressure on dust, for example, plays a major role in many models of the acceleration of winds from highly evolved, low-mass stars (e.g. MacGregor & Stencel 1992). The dust streams through the neutral gas and transmits momentum to it through collisions. The streaming of dust through neutral gas has also received attention in many other astrophysical contexts, including the evolution of dust-bounded H ii regions (Cochran & Ostriker 1977) and the radiation-driven implosion of dense globules (Sanford, Whitaker & Klein 1984). Radiation pressure on dust may levitate interstellar clouds above the disc of the Milky Way (Franco et al. 1991): in such clouds dust particles will stream through the neutral gas. No magnetic field was included in any of these studies. Hartquist & Havnes (1994) identified conditions under which dust grains are well-coupled to the magnetic field when the grains are driven by radiation pressure. In many circumstances the dust, gas phase ions, and electrons may be treated as a single fluid. Except for studies of the Wardle instability of shocks in dusty media (Wardle 1990; Stone 1997; Mac Low & Smith 1997), investigations of instabilities in weakly ionized astrophysical media on length-scales short compared to the Jeans length driven by fluids streaming relative to one another have been more limited.

There is a considerable literature on large-scale instabilities driven by the self-gravity of multifluid media (e.g. Mouschovias 1976; Nakano 1976; Huba 1990; Brandenburg & Zweibel 1995; Balsara 1996; Zweibel 1998; Kamaya & Nishi 2000; Mamun & Shukla 2001). These papers differ in aspects such as the number of different flow components assumed, the precise nature of the inter-species coupling and the inclusion of processes such as the self-gravity of the flow and large-scale gradients in flow properties. Many of these papers recover a large scale instability, first described by Mouschovias (1976) and Nakano (1976), in which the diffusion of magnetic field out of a self-gravitating clump reduces magnetic support, leading eventually to collapse. In some, rapidly growing, small length-scale instabilities are found (Huba 1990; Kamaya & Nishi 2000; Mamun & Shukla 2001), but, to date, rather restricted classes of relative orientation of magnetic field, mean flow and wavevector have been assumed.

We also note that there are many non-astronomical examples of interspersed multiphase flows, such as clouds, fluidized beds and microbial suspensions, which have been studied extensively. For example, Childress & Spiegel (1975) find buoyant instabilities, similar to those we discuss here, in systems of finite extent in both astrophysical and terrestrial contexts.

The past work on the Wardle instability and molecular cloud support has treated inhomogeneous media, in which the streaming is induced, for instance, by impulsive acceleration or by large-scale variations of the magnetic field. While astrophysical flows are necessarily inhomogeneous in the large, these variations can serve to obscure the mechanisms of small-scale instability. Given this wide variety of physical mechanism and equilibrium structure, in this paper we outline a general analysis for a simplified physical model, in which a charged magnetized fluid streams through a neutral fluid as a result of differential body forces. This model might most directly be related to flows with differential radiative forces on the fluids, but can be applied more widely. By assuming uniform initial conditions and treating the modes which we find as distributions, we can study the stability of short-wavelength modes in general, without needing to treat the specific global features which are important for longer wavelengths. Our analysis complements the previous work described above, by giving stability criteria for wave-vectors of arbitrary orientation and all initial angles between the body forces and magnetic fields, albeit for rather simpler input physics.

In Section 2, we present the basic two-fluid equations and derive the dispersion relation for linear waves. In Section 3, we present numerical solutions of the dispersion relation. For small wavelengths, we find that ‘resonances’ (where distinct modes have similar phase velocities) are important in understanding the stability properties, and discuss a graphical method of locating these resonances in general geometries. We then, in Section 4, analyse the stability of the solutions of the dispersion relation, proceeding from general analysis to specific analytic stability criteria for short and long wavelengths. These criteria compare well with the numerical results in the previous section, and confirm their generality. In Section 5 we apply the long-wavelength results to the properties of SiO maser spots in late-type stars. Finally, in Section 6, we summarize our results.

2 Basic equations and dispersion relation

In the present paper, we study the stability of two-fluid flows in which one component is coupled to a magnetic field. Differential forces on the two fluids lead to inter-phase slip in the equilibrium solution. We attempt to characterize the general properties which allow instabilities to feed off the slip energy. While the system we consider is simplified, it allows us to analyse the processes from which instability results in some detail.

The equations we treat are those of continuity and momentum for neutrals:  

1

 

2

and the ideal MHD equations for ions:  

3

 

4

 

5

where equation (5) maintains the solenoidal condition ∇·B = 0 so long as it is true initially, and the drag terms correspond to Stokes' law with drag coefficient λ. Here ρ_n, u_n, c_n and g are, respectively, the mean density, velocity, the effective isothermal sound speed and the net acceleration for the neutrals; ρ, u, c_i and g_i are the equivalents for the ions; and λ is a frictional coupling constant. The value of the vacuum permeability μ determines the magnetic field unit: conventional values include 1, 4π and 4π× 10⁻⁷. These equations are similar to those used in previous work (e.g. Shu 1983; Kamaya & Nishi 2000), except that in the present paper we do not restrict the relative orientation of the various vector fields (while, for the present, neglecting some of the physical terms included by these earlier authors).

Two distinct pressure terms are used for the two distinct phases. This might be taken as an assumption that the scattering between gas particles of the same phase is far more rapid than that between particles of differing phases. The limits c_i→ 0 and c_n→ 0 are relevant in particular contexts, but have degenerate eigenmodes: by assuming finite values of these parameters, the degeneracies are lifted within the present analysis.

The forces on the ionized component are often dominated by the effects of the overall curvature of the magnetic field, in which case one can assume u₀·B≃ 0 (Shu 1983; Mouschovias 1987). It is consistent to study small-scale instabilities in the present of such large scale gradients (e.g. Huba 1990). However, in this case at least part of the equilibrium force on one phase will not be proportional to mass, as we have assumed above: the detailed instability criteria will be somewhat different from those which we derive below, but the general behaviour should be similar. In what follows, we study the equations for general orientations of the various vector parameters, while bearing in mind the practical importance of cases in which u₀·B is small.

As in most previous papers, we have not included internal viscous terms for the individual phases, although collisional processes will in general lead to viscosity within the phases as well as inter-phase drag. This viscosity will lead to the stabilization of unstable wave modes at short wavelengths. We consider the effects of viscous terms in the context of an astrophysical example in Section 5, and verify that, in that case, viscosity can be neglected at the wavelengths which interest us. We will present a detailed treatment of the stability of flows including both inter-phase drag and internal viscosity in a future paper.

We perturb about an equilibrium with a constant slip velocity u₀ = u−u_n between the phases. The slip velocity satisfies the equilibrium condition  

6

Linearizing about this equilibrium with u_n = v_n, u = u₀ + v, ρ = ρ₀ + θ, ρ_n = ρ_{n, 0} + θ_n and B = B₀ + β, we find  

7

 

8

 

9

 

10

 

11

In what follows we will suppress indices 0 on ρ₀ and ρ_n,0. Looking for solutions of form exp i (k·x−ωt), we find  

12

 

13

 

14

 

15

 

16

These equations may now be manipulated to give a dispersion relation in the form , either directly or by noticing that they take the form of an eigenequation for eigenvalue ω, with eigenvector U^T = (θ_n, v^T_n, θ v^T, β^T).

Some factors of the dispersion relation can easily be found from equations (12)–(16). Taking the scalar product of k with (16), we find  

17

corresponding to an eigenvalue ω = u₀·k with eigenvector components β parallel to k and θ_n = v_n = θ = v = 0. This trivial eigenvector results from the requirement that ∇·B = 0 be maintained, and has an amplitude identically zero.

A further factor can be found, for which θ_n = θ = 0. Using equations (12)–(16), we find that  

18

with the components of the eigenvector having v parallel to B₀×k,  

19

 

20

In these modes, two Alfvén waves in the ionized gas are coupled to a shear mode in the neutral gas. We will refer to these as group A modes.

The remaining seven roots of the dispersion relation, which we will refer to as group B modes, are coupled by the remaining 2 components of the drag force. To simplify the form of the dispersion relation somewhat, we assume (without loss of generality) that the wavevector k is parallel to , and introduce scaled variables and parameters as follows:  

21

and a vector, v_A, with magnitude equal to the ionic Alfvén speed and direction parallel to B₀,  

22

Note that each of these variables has dimensions of velocity. This proves useful when interpreting the results of our analysis.

In terms of these variables, the dispersion relation for the group A modes, equation (18), takes the form  

23

The dispersion relation for the group B modes, obtained from equations (12)–(16), has the form

 

24

It is useful to write this latter dispersion relation in the form  

25

where  

26

 

27

 

28

 

29

 

30

 

31

and  

32

are the fast and slow ionic magnetosound wave speeds.

3 Numerical results

In this section, we present numerical results for the roots of the dispersion relations (23) and (25), concentrating in particular on circumstances where these roots correspond to physical instabilities. The roots were calculated as the eigenvalues of the complex matrix corresponding to equations (12)–(16), using the routine zgeevx from lapack (Anderson et al. 1999), and were verified by comparison with direct solutions to the polynomial dispersion relations. We show typical results in Figs 1 and 2.

Figure 1.

Variation of eigenvalues of group B with u_x, for various k. In this example, and u = (1, 0, −1) u_x. The axis Im(ω) = 0 is included to show where the solutions are unstable.

Figure 2.

Variation of eigenvalues of group B with k, for various u_x. As in Fig. 1, . The velocities u_x are chosen to select various mode resonances. The circled points show the short-wavelength limiting values. The axis Im(ω) = 0 is included to show where the solutions are unstable.

Where the coupling is weak (i.e. at high wavenumbers) and when the roots of the frozen system are well spaced, we find that the neutral acoustic modes are damped as ω=ω⁰− i λρ/2, the neutral shear modes as ω=ω⁰− i λρ and the ionized fast, slow and Alfvén modes all as ω=ω⁰− i λρ_n/2, in agreement with the results in Section 4.2.

Close to resonances, where two or more modes have similar phase velocities, additional modes need to be taken into account. The group A modes never become unstable, as confirmed by the analysis of Section 4.3. Instabilities are, however, found for group B modes. In Fig. 1, we show the values of for the group B modes, for one particular set of values of λ, ρ, c_i, ρ_n, c_n and B, and direction . The plots in this figure show the variation of the roots as the slip velocity is varied, while the value of the wavenumber k is increased through the set of plots. Numerous graphs have been included to illustrate the full range of topologies in phase space, and to allow the development of the instabilities to be followed from linear through non-linear order in λ.

For large k (or short wavelength), instabilities [where Im (ω) > 0] appear at resonances between the neutral sound waves and the ionized fast- and slow modes, and also (with a rather smaller growth rate) at resonances between one neutral shear mode and the ionized slow modes. As k decreases, the morphology of the phase space evolves, and additional modes begin to influence the instabilities of the resonant system. In particular, the roots corresponding to the slow modes of the undamped system merge at k∼ 25, leading to a cusped structure at Re (ω) = 0 for smaller k. The apex of this cusp corresponds to solutions at u_x→±∞. This structure is particularly notable since, as k decreases, the roots at increasingly large u_x become unstable (cf. Section 4.7). Nevertheless, comparing the panels of Fig. 1, it is apparent that the unstable modes which are seen can all be related to resonant instabilities for large k.

In a physical flow, waves with the full range wavenumbers k vary independently (at least in the linear limit), for fixed λ and u. In Fig. 2, we show the variation of the roots for varying magnitude of the wavevector, k, at fixed u_x (i.e. each of the points shown on these graphs appear simultaneously in a single physical system, as each corresponds to an independently-varying wave mode). As k decreases, each mode moves from its value at k→∞ (shown circled), at which its phase velocity is determined by the uncoupled equations. As these examples were all chosen with two wave modes close to resonance, the resonant modes show a characteristic behaviour for large values of k, diverging from their positions at k→∞ in a direction essentially parallel to the imaginary axis. In each case, as k→ 0 two of the solutions have Im (ω) →λ(ρ_i+ρ_n) = 11, while the remaining five solutions have Im (ω) → 0. For intermediate k, at least one root is on occasion unstable for each of the cases shown in Fig. 2(b)–(g). While each of these cases corresponds to a different mode resonance, and it is clear that these resonances have an important effect on the development of the roots with decreasing k, only three of the resonances lead directly to instabilities.

The properties of the most unstable modes in each of these cases are given in Table 1: the cases in which the resonant modes lead directly to instability have their most unstable modes at the largest wavenumbers. That the other cases are unstable is a result of the broad range of velocities away from resonance at which the three unstable mode couplings lead to instability, rather than resulting from the closer resonances which are present in these examples. Analytic criteria corresponding to these results will be determined in Section 4.4.

Table 1.

Properties of the most unstable solutions in Fig. 2. The panels of this figure corresponding to the various values of u_x are noted. Values are bracketed where the most unstable modes are not continuous with the resonant modes.

We have shown results for each of the values of k which are present for particular physical parameters, for a particular direction of the wavevector. In Section 4.1.2, we present a geometrical argument that allows us to generalize these results for the full range of wavevector directions present in a physical flow.

4 Stability analysis

In this section, we will study the stability of the solutions to the dispersion relations (23) and (25) analytically. First, we develop general stability criteria for dispersion relations with the general form of equation (25), including in particular a discussion of the case in which two or more of the undamped wave modes have similar frequencies. We then apply these results, and find that isolated modes are stable, and that the group A modes remain so for all λ. For the group B modes, however, we find that several two-mode resonances lead to instability. We also study higher resonances, and the stability of modes in the long-wavelength limit.

Our results confirm the importance of mode resonances for flow stability, as observed in the numerical results of the previous section.

4.1 General considerations

4.1.1 Perturbation theory

We first consider the limit in which the coupling parameter λ is small, i.e. that the damping parameters Ω_i and Ω_n in the dispersion relations are small. It is easy to see that P₀(Ω) may be written  

33

where Ω⁽⁰⁾_α = V_α are eigenvalues of the uncoupled system (i.e. when λ= 0, and so Ω_i = Ω_n = 0). We now use perturbation analysis, writing the frequency of each mode of equation (25) as a power-series expansion in the parameter λ, which is assumed small:  

34

where typically Ω⁽ⁿ⁾_α∝λⁿ.

If there are no degenerate or near-degenerate roots in the uncoupled system, then the first-order corrections may be written  

35

From this equation, it is apparent that for the flow to be stable, the roots of P₀(Ω) and Ω_iP₁(Ω) +Ω_nP₂(Ω) must interleave (cf. Appendix A; Whitham 1974). Since the roots of P₀(Ω) interleave with those of P₁(Ω) and P₂(Ω) individually, it is clear that, as Ω_i, Ω_n>0, this is indeed the case. [Note that this also demonstrates that Ω_iP₁(Ω) +Ω_nP₂(Ω) has six real roots, which will be important in our discussion of long-wavelength modes below.]

However, if there are two strictly degenerate modes of the uncoupled system, with Ω=Ω⁽⁰⁾_α, then equation (35) can no longer be applied and we have to take into account terms quadratic in λ in our perturbation analysis. In this case, equations (26)–(31) imply  

36

We then obtain the following quadratic equation for the first order correction to the degenerate mode, Ω⁽¹⁾_α,  

37

where P′₁(Ω) = d P₁/d Ω, etc. This equation has two solutions, which corresponds to the breaking of the degeneracy of the modes in the presence of a small perturbation. Note that even though the expansion was to second order in quantities proportional to λ, the correction is of first order.

If we have near-degenerate roots, then V_α∼V_β for some set of roots {α, β, …}: we will refer to this circumstance as a resonance. As a result, we have to consider the difference(s) V_α−V_β, α≠β as further small parameter(s). For a two-wave resonance, we have  

38

where the root of the quadratic with the smaller real part should be chosen for continuity with the non-resonant case.

This equation can be written in the form  

39

with  

40

and  

41

We can now apply the Hermite-Biehler theorem (see Appendix A), which tells us that the system described by the equation  

42

is stable (that is the roots of the equation have negative imaginary parts) when the roots of the polynomials Q_n and Q_n−1 are real and interleave (and their leading terms are positive), and that if this is not the case the system is unstable. In the case of equation (39) we have polynomials of first and second degree, and can easily obtain conditions for stability. For higher resonances, the expressions are clumsy, so it is better to find instability conditions for particular cases rather than apply the general expression.

4.1.2 Geometry of resonances

In the preceding section, we presented a set of general criteria for the stability of flows with internal damping parameters, and suggested that mode resonances may play an important role. We now find general conditions for the presence of mode resonances in two-fluid MHD flows in three dimensions, by comparing the geometry of phase diagrams for the magnetosound modes of the ionized gas component and the sound and shear modes of the neutral component.

In Fig. 3(a) and (b), we show the phase speeds for the fast, slow and Alfvén modes, as a polar plot of ω/k as a function of the angle between k and the magnetic field B, while in Fig. 3(c) and (d) we show the same for sound waves and a shear mode in a frame in relative motion as a function of the angle between k and u₀.

Figure 3.

Phase velocity profiles, v =ω/ k, as a function of the angle between k and B or u₀, plotted on polar axes. (a), (b) Solid — fast and slow, dashed — Alfvén, [the dotted curve — sound — is included for comparison but has no dynamical role here], c_i > v_A for (a) and c_i < v_A for (b); (c), (d) Solid — sound, dotted — shear, with velocity offset, (c) is for subsonic flow and (d) for supersonic flow.

It is clear that the passing of the phase curve for sound waves in Fig. 3(d) through the origin corresponds to the Cerenkov condition for the emission of sound waves by an individual particle moving at u₀. In a similar manner, resonance conditions between modes in the different phases can identified by looking for intersections between curves on suitably scaled and oriented versions of these plots (remembering to keep the origin at the same place in the ionized and the neutral plot). For example, resonances will always occur between the slow and shear modes, unless u₀ is parallel to B and u₀ > v_A. Likewise, resonances between the slow and neutral sound modes will be present in general once the relative motion is faster than the neutral sound speed, and for some choices of parameters once the slip velocity is greater than the difference between the neutral sound speed and the Alfvén velocity.

For finite damping, the resonance condition is weakened: the lines shown in Fig. 3 may be thought of as blurred, although as a side-effect additional modes must also be considered in the stability analysis.

We have now determined general stability criteria, suggested the importance of mode resonances in determining stability, and demonstrated that they will be an almost universal feature of MHD flows with drift velocities between components. However, depending on the form of the frictional force, only some of these resonances will result in flow instabilities. In the following subsections, we will apply these results to the dispersion relations given in Section 2 to determine quantitative stability criteria.

4.2 Non-resonant modes

Using relation (35), first let us consider the stability of the modes of equation (25) when there are no resonances in the uncoupled system, and Ω_i and Ω_n are small. We find that all the modes are stable. Details are given below.

Neutral shear modes For neutral shear modes, with Ω⁽⁰⁾_α = 0,  

43

Neutral sound waves For neutral sound waves, with Ω⁽⁰⁾_α = ±c_n,  

44

Magnetosound and Alfvén modes for the fast and slow magnetosound waves, with Ω_α⁽⁰⁾ = u_x±v_f,s , and also for the (group A) Alfvén modes, with Ω⁽⁰⁾_α = u_x±v_Ax ,  

45

4.3 Group A modes

The stability of the group A solutions can easily be analysed for all λ (not necessarily small), by applying the Hermite-Biehler theorem to the dispersion relation, equation (23). This equation couples two Alfvén modes and one neutral shear mode, and can be written  

46

where  

47

and  

48

It is clear that both Q₂ and Q₃ have their full complement of zeroes, and that the zeroes of Q₂ interleave with those of Q₃ (so long as there are no strict degeneracies and Ω_i, Ω_n>0). Hence, by the Hermite-Biehler theorem, the group A modes are always stable.

4.4 Two-wave resonances

We now consider the stability of the flow close to resonances between various pairs of group B modes. In order to determine whether instability occurs, we will follow the method of Section 4.1.1.

4.4.1 Ionic magnetosound—neutral sound resonance

First, we consider the case of a resonance between a neutral sound mode and an ionic magnetosound mode, so  

49

where α and γ can independently take values +1 or −1. Applying the Hermite-Biehler theorem to the corresponding equations (39)–(41), one finds that instability occurs when the following inequalities hold:  

50

and  

51

In the limit of strict degeneracy, the right-hand side of the second condition goes to zero, so in this limit the condition is no constraint. Note also that equation (50) is the condition that the right-hand side of equation (51) is non-negative. Equation (51) is in effect a minimum condition on λ/k for instability to occur (implying instability for large damping or long wavelength), for modes at a given distance from resonance. Note that these properties of conditions (50) and (51) also apply to conditions (57) and (58) below.

We can now consider particular cases of these general relations.

• (i)

  αc_n∼u_x+γv_f.

  So long as v_Ax≠ 0, equation (50) may be written  

  

  52

  If v_A⊥u then  

  

  53

  If c_i = 0 as well, then  

  

  54

  which is possible only if γα=−1, i.e. when the sound and fast-mode waves are oppositely directed.
• (ii)

  αc_n∼u_x+γv_s.

  Equation (50) gives us  

  

  55

  If v_A ⊥ u then  

  

  56

  which is possible only if γα= 1 , i.e. where the phase velocities of the sound and slow-mode waves are in the same direction. Unlike the previous case, we cannot take c_i = 0 in this relation, because in that case the remaining slow-mode wave will also have a similar frequency. The resulting higher resonance will be considered below.

4.4.2 Ionic magnetosound—neutral shear mode resonance

We now consider the case when u_x+γv_f,s∼ 0. Applying the Hermite-Biehler theorem to equations (39)–(41) gives the instability conditions  

57

and  

58

As before, we consider two particular cases

• (i)

  u_x +γ v_f ∼ 0.

  In this case equation (57) gives us  

  

  59

  If v_A⊥u then the interacting modes are always stable.

• (ii)

  u_x +γ v_s ∼ 0.

Now we have from (57) 

60

so that if v_A⊥u then this mode interaction will always lead to instability for sufficiently large λ, so long as no other modes start to interact.

4.5 Degenerate slow magnetosound waves

If we have c_i∼ 0 (which is often the case in astrophysical applications), then the two slow magnetosound waves are resonant, i.e. u_x−v_s∼u_x+v_s because v_s∼ 0. This resonance is of a different nature to those discussed in Section 4.4, as it is between two modes that propagate in the same phase.

When c_i is strictly zero, two slow magnetosound waves become degenerate with common eigenvalue Ω⁽⁰⁾_α = u_x. In order to investigate stability in this case, we have to solve equation (37), and find  

61

Instability can occur for the mode with positive choice of sign if  

62

For small but finite c_i, for instability to occur then in addition to condition (62) we also require that  

63

For the important case v_A⊥u, condition (62) leads to the simple relation  

64

so long as v_Ax≠ 0. Note that the perturbation analysis will not apply when u²_x is close to c²_n, since Ω⁽¹⁾_α given by equation (61) will be no longer be small. While no resonance with modes in the neutral phase has been assumed, this condition clearly relates to a resonance effect with the neutral sound waves in some fashion (indeed, when followed in numerical solutions, the mode is continuous with instabilities which originate at the slow magnetosound/neutral sound resonance for finite c_i). As the ionized sound speed becomes small, the system of coupled slow magnetosound modes becomes sensitive to the presence of other modes, even those rather far from apparent resonance conditions (cf. the sensitivity of the slow-mode waves to driving by other modes of non-linear amplitude in the ionized gas, Falle & Hartquist 2002).

4.6 Higher resonances

We now consider the stability of some higher order mode couplings.

4.6.1 Ionic slow magnetosound-neutral shear mode resonance

Let us consider the case when u_x−v_s∼ 0 ∼u_x+v_s, which is possible if c_i and u_x are small.

First let us investigate the case c_i = u_x = 0, in which when we have a three times degenerate root Ω⁽⁰⁾_α = 0 for the uncoupled system. If we apply perturbation analysis to these roots, we obtain the following equation for the correction Ω⁽¹⁾_α:  

65

The polynomial Q₂(Ω⁽¹⁾_α) has two identical roots Ω⁽¹⁾_α = 0, which means that in order to satisfy the interleaving criterion, the polynomial Q₃(Ω⁽¹⁾_α) must also have one root Ω⁽¹⁾_α = 0. This is possible only if v_A·u = 0– indeed, for  

66

Q₃ has only a single real root. From this we conclude that the system is marginally stable when v_A·u = 0 and unstable when v_A·u≠ 0. In the latter case, the equation has one term of second order, so the roots are  

67

which scale as λ^2/3 for λ→ 0, with one root being unstable.

Now we consider the general resonant case where u_x and c_i are small but finite. For the Ω⁽⁰⁾_α = 0 mode we have  

68

The corrections derived for the roots which have Ω⁽⁰⁾_α = u_x+γv_s give equivalent results, relative to their alternative datum, as must be the case since the three roots of equation (68) correspond to the three resonant modes.

The Q₂ polynomial for equation (68) has two real roots  

69

For the Q₃ polynomial, we first treat the case v_A·u = 0. Here, there are three real roots, one of which is zero and other two are  

70

Interleaving of these roots with those given in equation (69) is violated (that is, the system is unstable) when  

71

so we see that the system, which is marginally stable under fully degenerate conditions, can be unstable under resonant conditions. The range of unstable conditions widens in velocity as Ω_i and Ω_n increase.

If v_A·u≠ 0, we have an equation of third degree for Q₃,  

72

Again the last term in this equation is of second order, so in the short-wavelength limit the roots are still given by equation (67).

In the more general case where Ω_iΩ_n ≳ v²_A(v_A·u)²/v⁴_Ax, equation (72) has complex roots (so the system is unstable) if  

73

(from the usual relations for cubic equations, e.g. Press et al. 1992), where we have introduced new variables  

74

 

75

This implies that the solutions are unstable as v_A·u→±∞, other terms being equal: in particular, if u_x is zero there are complex roots when  

76

If ∣Γ₂ȣ≪Γ^3/2₁, condition (73) is satisfied for u_x in the range  

77

If inequality (73) is violated, then equation (72) has three real roots and we have to solve the problem of interleaving. Criteria for interleaving may be obtained by inserting the roots of the Q₂ polynomial, equation (69) into the left-hand side of equation (72). These are rather complex criteria to apply, except numerically. In general, however, for small Ω_i, Ω_n the flow is unstable for conditions along whichever of the lines u_x = ±v_s satisfies equation (60); as Ω_i, Ω_n increase this region of instability broadens and also expands along the other of u_x = ±v_s, in agreement with the limiting cases given by equations (71) and (77).

4.6.2 Neutral sound—neutral shear mode resonance

Let us consider the case αc_n∼ 0. One can show that equation for the corrections has the form  

78

It is obvious that roots of Q₃ and Q₂ always interleave, and so the system is always stable under these resonant conditions (as is clear physically).

4.6.3 Five-wave resonance

There is also the possibility of a five-wave resonance if 0 ∼αc_n∼u_x+γv_s. However this case gives us a fifth-order polynomial, which would be extremely complicated to analyse. We therefore restrict ourselves here to investigating only the case of full degeneracy, e.g. c_n = 0, c_i = 0 and u_x = 0. Then we have three roots Ω= 0, and for the remaining two we have the equation for the corrections,  

79

It is clear from this that the system is stable. However, as we have seen above, the marginal stability of the three roots Ω= 0 under degenerate conditions does not guarantee stability under resonant conditions.

4.7 Long-wavelength stability

We now consider the stability of modes in long wave limit, i.e. for small wavenumbers. Note first that if we take the limit k→ 0 in equation (24) then five modes have behaviour ω→O(k) and two modes have behaviour ω→−i λ (ρ_i+ρ_n) +O(k)[cf. also the behaviour of the roots for long wavelengths shown in Figs 1(l) and 2].

Now we have to find the corrections to these limits of order k:  

80

and  

81

According to the definition of the limit, small k means that the second term in equation (81) is much smaller then the first one. As the first term has a negative imaginary part, the smaller second term cannot change the overall sign. Hence the roots given by equation (81) are stable in the long-wavelength limit.

It is more difficult to determine the stability of the roots given by equation (80). To first order in k, the stability of these five modes in the limit of long wavelengths is determined by the roots of  

82

together with the condition that the leading term in this equation has a positive coefficient (which is the case). This equation is independent of k, and its roots are the velocities V_α in equation (80). As equation (82) is a polynomial equation with real coefficients, complex roots can arise as conjugate pairs and any such complex roots in equation (82) will result in instability in the limit of small k.

Real roots of equation (82) are marginally stable to first order. To determine their asymptotic stability, we would need to study the interleaving of these roots with those of Ω_iP₁(Ω) +Ω_nP₂(Ω): if the roots of these two polynomials do not interleave, the asymptotic solution will be unstable at order k². Note that we know that Ω_iP₁(Ω) +Ω_nP₂(Ω) has its full complement of real roots, from the discussion after equation (35) above.

We will first study the stability of the long-wavelength solutions in various limiting cases. For u = 0 we have five real (and hence stable) roots, which correspond to the fast- and slow-mode waves of the fully coupled system together with an additional mode with Ω= 0, as would be expected on physical grounds.

If u_x is large compared to any other characteristic velocity, the solutions are Ω= 0[or more accurately Ω_nc²_n/(Ω_iu_x)], and Ω_iu_x/(Ω_i+Ω_n) and u_x (the latter two are double roots). This can be easily seen if we write equation (82) as follows  

83

The presence of double roots means these solutions are only marginally stable in the limit, so we need to carry the analysis to higher order to determine stability in the asymptotic regime. To lift these degeneracies, we should expand to lower order in Ω and u_x, or to higher order in k (or include non-linear terms neglected in our initial linearization of the problem).

If k⊥v_A, i.e. v_Ax = 0 with all other terms finite, equation (83) has roots Ω = u_x and Ω = Ω_iu_x/(Ω_i+Ω_n). The remaining three roots of the equation satisfy  

84

where α_i = Ω_i/(Ω_i+Ω_n), α_n = Ω_n/(Ω_i+Ω_n). Considering the ordering of the roots of the factor proportional to α_i and that proportional to α_n≡ 1 −α_i, it is clear that there must be at least two real roots of equation (84) for all 0 < α_i < 1, and hence there will be three real roots. Therefore for v_Ax = 0, all the solutions are stable to order O(k).

If we consider corrections to the O(u_x) asymptotic solutions, then for the two roots Ω = u_x+δ, we find that the highest order corrections have δ = O(1/u_x), and are given by  

85

which can clearly be complex. For example, if v_A·u = 0 and c_i = 0 the criterion for instability becomes 4v²_Axc²_n/v⁴_A>Ω_i/Ω_n. Similar expressions for the case Ω≃Ω_iu_x/(Ω_i+Ω_n) are more complex, but the condition for instability is  

86

If, instead, we include higher order terms in k, where in the limit of large u_x 

87

 

88

 

89

then the common roots Ω= 0 and u_x (twice) are marginally stable. Interleaving is satisfied for the remaining four roots, so they are each stable for all k at large u_x (although only the limiting form as k→ 0 is required for the present discussion).

First-order instabilities can also arise in the intermediate-velocity regime. Looking again at equation (83), we have shown that quintic consisting of the first two terms has five real roots. However, the values v_Ax and v_A·u can be varied independently of any of the other parameters of the overall equation, so the additional linear function consisting of the third and fourth terms is entirely arbitrary. Even when v_A·u = 0, the number of real roots of equation (83) can be as few as one.

As the long wave limit is quite important for the astrophysics of molecular clouds and stellar winds, in the next section we will study a typical astrophysical example, SiO maser spots in late-type stellar winds.

5 SiO maser spots in late-type stellar winds

As mentioned in the introduction, the current results are applicable to mass loss from highly evolved stars, which are believed to be driven by radiation pressure on dust particles (Bowen 1988; MacGregor & Stencel 1992; Mastrodemos, Morris & Castor 1996; Simis, Icke & Dominik 2001). The winds are subject to strong perturbations as a result of stellar pulsation, and models have found both radiation-pressure instabilities and others resulting from the condensation of the dust grains, as well as radiative instabilities. The structures which result from these instabilities may be important in determining the properties of planetary nebula halos which will form as the star ages further. Mastrodemos et al. (1996) find that slip between dust and gas components leads to formation of dense shells at intermediate radii, but that these dissipate as the wind moves away from the star. Simis et al. (2001), however, find that dense shells survive to large radii, and that an accurate treatment of dust-gas slip is essential in treatments of late-type winds.

If, in the wind of a late-type star, ion-neutral friction is the main source of coupling and the streaming speed is less than about 30 km s⁻¹, λ can be calculated from the results obtained by Osterbrock (1961) for the ion-neutral momentum transfer cross-section, so  

90

where x is the fractional ionization, i.e. the ratio of the number density of ions to the number density of hydrogen nuclei. The value of k in the dimensionless units used in Figs 1 and 2 is related to the wavenumber in physical units k_phys by k = k_phys/k₀, where  

91

To compare this with the scalelengths of features in the winds, we must have some idea of the appropriate values of v_Ai, n_n, x and of the scalelengths of interest. SiO masers are found in the outflows of some evolved stars at positions near those at which dust is expected to form (e.g. Kemball & Diamond 1997). The magnetic field strength in the maser spots can be inferred with difficulty from the observed polarization (e.g. Watson & Wyld 2001), although it is not certain that the polarization is in fact caused by the Zeeman effect (Wiebe & Watson 1998). The resulting estimates of the magnetic field strengths are in the range of 2 to 10 G (Elitzur 1996; Kemball & Diamond 1997). Pumping models for these masers suggest that the neutral density in the spots is n_n≈ 10¹⁰–10¹¹ cm⁻³ (Doel et al. 1995). If 1 per cent of the mass is contained in grains, v_Ai could be as high as 10³ km s⁻¹, but may be closer to 100 km s⁻¹. The size of individual maser spots is ≃10¹² cm (e.g. Kemball & Diamond 1997), for which k≈ 3 × 10⁻¹³ cm⁻¹ if the maser spot size corresponded to a half wavelength of a perturbation. The fractional ionization is uncertain, but it is unlikely that grains will form where carbon is ionized, and some of the elements with low ionization potentials will be depleted substantially from the gas phase due to the grain formation process. Grain—neutral friction for a grain—neutral relative speed comparable to the thermal speed of the neutral material is of similar magnitude to the ion-neutral friction for x between 10⁻⁷ and 10⁻⁶, if the grains have a fractional abundance and size distribution similar to those in the interstellar medium (e.g. Baker 1979).

Therefore the SiO maser regions are probably large enough that their overall properties correspond to our long-wavelength limit. To study their stability in detail, we consider parameters characteristic of the regions of late-type stellar winds with SiO maser spots, as follows: ρ_i/ρ_n = 0.01, c_n = 3 km s⁻¹, c_i = 0, v_A = 1000 km s⁻¹ and u = 0 − 100 km s⁻¹. We also assume that u is perpendicular to magnetic field, i.e. v_A·u = 0, and in particular study the case where the wavevector, k, bisects the angle between u and v_A. Numerical solutions of the long-wavelength dispersion relation are shown in Fig. 4. For comparison we also show in Fig. 5 the calculations for smaller Alfvén speed v_A = 20 km s⁻¹. Note that for some slip speeds there are two pairs of complex roots, but that for larger Alfvén velocities the higher velocity roots are stable at all slip speeds. We will refer to the instability of the higher velocity roots, shown in Figs 5(a) and (b), as a type I instability, and that in Figs 4, 5(c) and (d) as type II.

Figure 4.

Values of Ω for limit of long wavelengths, for v_A = 1000 km s ⁻¹. (a) and (b) show the real and imaginary parts of two roots which are complex for part of the velocity range. The remaining three roots are stable and not shown here. The roots do not tend to their asymptotic behaviour in this plot, as the maximal slip velocity u_max = 100 is small relative to the Alfvén velocity (cf. the case of v_A = 20 km s ⁻¹ in Fig. 5 ).

Figure 5.

Values of Ω in the long-wavelength limit, for v_A = 20 km s ⁻¹. (a) shows the real parts of the roots for small u_x. The regions in (a) where the real parts are distinct correspond to stable roots with zero imaginary parts: (b) shows the imaginary parts of the complex conjugate pairs. (c) and (d) are the same as (a) and (b), except that they are plotted for a rather wider range of u_x, and the asymptotic behaviours of the roots are shown as dashed curves. At large u_x, the roots with the largest real parts follow the expected asymptotic behaviour, but have finite imaginary parts. The next two roots should converge to a common asymptotic limit, Ω≃ 0.01u_x , but while the variation of the smaller of this pair is almost indistinguishable from this form (indeed, it overlies the dashed curve showing the limit on the graph), the larger lies well away from its asymptotic value, and convergence is only clear at far larger u_x , as shown in (e).

It will be noticed that as v_A changes from 20 to 1000 km s⁻¹, the region of type II instability moves to higher slip velocities. For v_A = 20 km s⁻¹, the maximum of Im (Ω) is around u_x∼ 1.5 km s⁻¹ and for v_A = 1000 km s⁻¹ this maximum is around u_x∼ 20–25 km s⁻¹. The value of the maximum growth rate, however, changes little, and is around 0.06–0.1 km s⁻¹.

For small Alfvén velocities (like 20 km s⁻¹), a type I instability is also present. The regions of two instabilities overlap to a degree, with the range of type II instability being for u_x≃ 1.2 − 1.6 km s⁻¹ and that for the type I instability being for u_x≳ 1.5 km s⁻¹.

The maser spots are typically at distances of 5–10 au from the centres of the stars. For an outflow speed of the order of 10 km s⁻¹, we conclude that an instability will have no significant effect on the maser spots unless it grows on a time-scale of a few years or less. From Fig. 5(b) one can see that for relatively small Alfvén velocities this is indeed the case, and such maser spots would be destroyed by type I instability. But as maser spots exist at these radii, we can conclude that magnetic field (and therefore the Alfvén velocity) cannot to be too small. Indeed, for Alfvén velocities around 1000 km s⁻¹, type-I instability is suppressed and we only have to deal with type-II instability (see Fig. 4). This instability at its maximum grows on time-scale of order 5–10 yr, but the range of streaming velocities at which the instability has even this large a growth rate is quite narrow, between 28 and 35 km s⁻¹. The streaming speed of grains through a maser spot is uncertain and depends on the fractional ionization and the grain size distribution; it may be that in reality the streaming speed is too large or too small for the instability to be driven at a rate high enough to destroy the maser spot. It will be important, as models of late-type stellar wind including both the chemistry of grain formation and dynamics are developed (e.g. MacGregor & Stencel 1992; Cherchneff 1998; Gail & Sedlmayr 1998), to determine the relationship between the calculated grain—neutral streaming speed and the input magnetic field. These results are necessary for a full application to stellar outflows of the present analysis.

Finally, it is of interest to determine whether the short-wavelength analysis developed in previous sections can be applied to the finer-scale properties of the maser spots. This will be the case if the widths of the resonant interactions are narrower than the spacing of non-resonant modes, i.e. that  

92

where V_min is the smallest characteristic velocity of the problem, which should be u∼ 1 km s⁻¹, and we take ρ_n in the left hand side as it is larger than ρ_i. This inequality means that, in this example, short wavelengths are 2π/k≲ 10¹⁰ cm.

This maximum wavelength must be compared to minimum wavelengths required for the fluid approximations to be valid, and for the growth of the two-fluid instability to outweigh wave damping as a result of viscosity. The maximum mean free paths for both electrons and ions are  

93

(Melrose 1986), while for neutral atomic species, the free path is  

94

 

95

assuming a geometric cross-section. For longer wavelengths, the finite collision frequencies will lead to viscous damping of waves in a number of periods given by the ratio between the wavelength and the free path, and any instability must grow rapidly enough to outweigh this damping. For a given slip velocity, the growth rates for the instability we discuss are proportional to the density while the viscous damping rates are proportional to the ratio of the square of the frequency divided by the density, so damping will indeed only be important for high-frequency waves in diffuse media. In fact, at these smaller scales further instability modes may also become important (driven by, e.g. the Hall effect, Balbus & Terquem 2001).

For the particular example under discussion, we require 10⁵ cm ≪λ_w≪ 10¹⁰ cm which obviously can be satisfied, so that our short-wavelength analysis is applicable for a substantial range of wavelengths. In general, an appreciable range of potentially unstable wavelengths remains for all temperatures T≪ 10⁶ K. While the shorter wavelengths have less direct observational significance, their growth to non-linear amplitudes will likely have substantial effects on the overall flow structure.

6 Discussion and conclusions

We have seen that the slip between different fluids in a multifluid medium can drive instabilities. At short wavelengths, several modes of instability result from resonances between waves propagating in the different fluids; these instabilities remain for longer wavelengths but they become increasingly intercoupled and less easy to characterize. In particular, neutral shear/ionized slow-mode resonance instabilities will grow in almost all cases where there are finite drift velocities. The maximum growth rates we find are of order the characteristic ion/neutral collision rates in gas without internal slip, which are typically  

96

(Osterbrock 1961; Brandenburg & Zweibel 1995), rather than inversely proportional to density as is the case for instabilities derived from leakage of magnetic flux as a result of ambipolar diffusion. This corresponds to a characteristic length-scale of 2.5 × 10⁻⁴(v/1 km ⁻¹)n⁻¹_n pc, where we note that the most unstable wavelengths may be one hundredth of this (Table 1). The instabilities appear to be similar to the two-stream instability of plasma flows (e.g. Melrose 1986; Bingham et al. 2000).

In Section 5, we apply our analysis in its long-wavelength limit to the overall properties of SiO maser spots in late-type stellar winds. We find that if the SiO spots are not to be subject to violent instabilities, the slip velocity between the phases in these regions must lie outside certain limits, or the magnetic field in these spots must be very strong.

Our short wavelengths results also have important implications for many other astrophysical systems. In molecular clouds, we expect that the ionized sound speed is small and that the magnetic field is perpendicular to the inter-component slip velocity, but that this slip velocity is not small. This means that we always can find directions such that inequality (71) is satisfied. It seems likely, therefore, that molecular clouds are generically unstable to the growth of slow-mode waves. In the present analysis we can only conjecture the non-linear endpoint of these instabilities (although Franqueira, Tagger & Gómez de Castro 2000 find a non-linear filamentation process in two fluid MHD turbulence which may be related to the instabilities we describe). The two obvious possibilities — fractionation of the phases as a ‘slugged’ flow, with consequent rapid loss of magnetic field support, or the limiting of the wave spectrum at finite amplitude — each have clear practical and observational consequences for the ecology of the interstellar medium.

Our present study is substantially simplified. This has allowed us to derive some rather general results. In future work, we will model the structure of the systems of interest more completely, including the spatial structure of the background flow, more interacting phases (e.g. treating electrons, ions, neutrals and a spectrum of sizes of dust particle as independent species), variation of the frictional constants with state, and the detailed non-linear evolution of the instabilities.

Acknowledgments

PVT thank the Royal Society and PPARC for support during this work. RJRW is supported by a PPARC Advanced Fellowship, and thanks the Department of Physics and Astronomy in Leeds for hospitality while this work was developed. We wish to thank Tom Hartquist for his role in suggesting this work, his continuing interest in its development and useful comments on the manuscript. Helpful comments from the anonymous referee and from J. Franco were also much appreciated.

Appendix

Appendix A: The Hermite-Biehler theorem

We quote here without proof the Hermite-Biehler theorem, a very general result which underlies much of our analysis. From theorem 4′ of Levin (1964), Chapter VII:

Theorem 1. (Hermite-Biehler) In order that the polynomial  

(A1)

where u(z) and v(z) are real polynomials, does not have any roots in the closed upper half-plane Im (z) ⩾ 0, it is necessary and sufficient that the following conditions are satisfied.

• (i)

  The polynomials u(z) and v(z) have only simple real roots, and these roots separate one another, i.e. between two successive roots of one of these polynomials there lies exactly one root of the other;

• (ii)
  at some point x₀ of the real axis,  

  

  A2

  where u′(z) = d u(z)/ d z and v′(z) = d v(z)/ d z.