MHD Waves in Weakly Ionised Media

Hydromagnetic Waves in Weakly Ionised Media. I. Basic Theory, and Application to Interstellar Molecular Clouds


We present a comprehensive study of MHD waves and instabilities in a weakly ionised system, such as an interstellar molecular cloud. We determine all the critical wavelengths of perturbations across which the sustainable wave modes can change radically (and so can their decay rates), and various instabilities are present or absent. Hence, these critical wavelengths are essential for understanding the effects of MHD waves (or turbulence) on the structure and evolution of molecular clouds. Depending on the angle of propagation relative to the zeroth-order magnetic field and the physical parameters of a model cloud, there are wavelength ranges in which no wave can be sustained as such. Yet, for other directions of propagation or different properties of a model cloud, there may always exist some wave mode(s) at all wavelengths (smaller than the size of the model cloud). For a typical model cloud, magnetically-driven ambipolar diffusion leads to removal of any support against gravity that most short-wavelength waves (or turbulence) may have had, and gravitationally-driven ambipolar diffusion sets in and leads to cloud fragmentation into stellar-size masses, as first suggested by Mouschovias more than three decades ago – a single-stage fragmentation theory of star formation, distinct from the then prevailing hierarchical fragmentation picture. The phase velocities, decay times, and eigenvectors (e.g., the densities and velocities of neutral particles and the plasma, and the three components of the magnetic field) are determined as functions of the wavelength of the disturbances in a mathematically transparent way and are explained physically. Comparison of the results with those of nonlinear analytical or numerical calculations is also presented where appropriate, excellent agreement is found, and confidence in the analytical, linear approach is gained to explore phenomena difficult to study through numerical simulations. Mode splitting (or bifurcation) and mode merging, which are impossible in single-fluid systems for linear perturbations (hence, the term “normal mode” and the principle of superposition), occur naturally in multifluid systems (as do transitions between wave modes without bifurcation) and have profound consequences in the evolution of such systems.

diffusion — ISM: magnetic fields — MHD — plasmas — stars: formation — waves

1 Introduction – Background

A typical molecular cloud which has not yet given birth to stars is a cold () but complex, partially ionised system, in which self-gravitational and magnetic forces are of comparable magnitude, with thermal-pressure forces becoming important at high densities () or along magnetic field lines. Mouschovias (1976) showed that, barring external disturbances, if the magnetic field were to be frozen in the matter, interstellar clouds that have not yet given birth to stars would remain in magnetohydrostatic (MHS) equilibrium states. However, ambipolar diffusion (the relative motion of neutral particles and charged particles attached to magnetic field lines) is an unavoidable process in partially ionised media. It reveals itself in two distinct ways, depending on whether it is magnetically or gravitationally driven (see discussion in § 4). The two kinds of ambipolar diffusion acting together initiate fragmentation and star formation in molecular clouds (Mouschovias 1987a).

In this fragmentation theory, the evolutionary (or fragmentation, or core formation) timescale is the gravitationally-driven ambipolar-diffusion timescale, . This does not mean, however, that it takes a time equal to to form stars. The star-formation timescale can be a fraction or a multiple of , depending on the mass-to-flux ratio of the parent cloud and the degree to which hydromagnetic waves (HM) contribute to the support of the cloud: the closer to its critical value the mass-to-flux ratio is and/or the greater the contribution of HM waves to cloud support, the faster the evolution and the shorter the star-formation timescale (e.g., see Mouschovias 1987a; Fiedler & Mouschovias 1993, Fig. 9a; Ciolek & Basu 2001; Tassis & Mouschovias 2004, Fig. 4). 1 The collapse retardation factor, (where is the free-fall time and the neutral-ion collision time), is the factor by which magnetic forces slow down the contraction relative to free-fall. This was a new theory of fragmentation (or core formation), initiated by the decay, due to ambipolar diffusion, of relatively small-wavelength perturbations (Mouschovias 1987a, 1991a).

Magnetic braking operates on a timescale shorter than the ambipolar-diffusion timescale and even the free-fall timescale, and keeps a cloud (or fragment) essentially corotating with the background up to densities and thus resolves the angular momentum problem of star formation. More specifically, the entire range of periods of binary stars from 10 hr to 100 yr was shown to be accounted for by this self-initiated mode of star formation (Mouschovias 1977). Even single stars and planetary systems become dynamically possible (Mouschovias 1978, 1983).

Star formation, whether self-initiated or triggerred (e.g., by a spiral density wave, or the expansion of an Hii region or a supernova remnant; see review by Woodward 1978) is an inherently nonlinear process. Ambipolar-diffusioninitiated star formation has been studied analytically (Mouschovias 1979, 1991a, b) and numerically using adaptive grid techniques in axisymmetric geometry up to densities , by which isothermality begins to break down (Fiedler & Mouschovias 1992, 1993; Ciolek & Mouschovias 1993, 1994, 1995; Basu & Mouschovias 1994, 1995a, b). More recently, these calculations were extended into the opaque phases of star formation (Desch & Mouschovias 2001; Tassis & Mouschovias 2007a, b, c; Kunz & Mouschovias 2009, 2010). The key conclusions of the earlier analytical calculations have been verified and numerous new, specific, quantitative predictions have been made, many of which have been confirmed by observations (e.g., see Crutcher et al. 1994; Ciolek & Basu 2000; Chiang et al. 2008; reviews by Mouschovias 1995, 1996). One may nevertheless take a step back from the relatively complicated numerical calculations and ask which results, if any, of the nonlinear simulations can be recovered with a linear analysis. With one’s confidence increased in the validity of the linear approach (within some self-evident limits), one may then make new predictions concerning phenomena that cannot be or have not been included yet in the nonlinear calculations.

The propagation, dissipation, and growth of perturbations in a physical system depends both on the nature of the perturbations and the properties of the system. Hydromagnetic waves (or MHD turbulence) seem to play a significant role in molecular clouds on lengthscales typically greater than . They have been shown to account quantitatively for the observed supersonic but subAlfvénic spectral linewidths: an observational almost-scatter diagram of linewidth versus size is converted into an almost perfect straight line if plotted in accordance with a theoretical prediction by Mouschovias (1987a), which relates the linewidth, the size, and the magnetic field strength of an observed object (see Mouschovias & Psaltis 1995, Figs. 1 and 2, and update by Mouschovias et al. 2006, Figs. 1 and 2).

Using a linear analysis as a first step in understanding nonlinear phenomena is not, of course, a new idea. Jeans (1928) used it to obtain his famous instability criterion for the collapse of a cloud against thermal-pressure forces. Hardly any astrophysical system exists whose stability with respect to small-amplitude disturbances has not yet been studied by using at least an idealized, mathematically tractable model of the physical system. A magnetically supported molecular cloud, however, defies a simple linear analysis. First, no realistic equilibrium states have been obtained by analytical means. Second, to study the role of ambipolar diffusion in star formation, one must use at least the two-fluid magnetohydrodynamic (MHD) equations governing the motions of the neutral particles and the plasma (ions and electrons). In fact, as shown by Ciolek & Mouschovias (1993), charged (and even neutral) grains play a very significant role in the ambipolar-diffusioninitiated protostar formation. One then has to use at least the four-fluid (neutral molecules, plasma, negatively-charged and neutral grains) MHD equations even for a linear analysis to be realistic and relevant to typical molecular clouds. In this paper we use the two-fluid MHD equations to study the propagation, dissipation, and growth of HM waves in an idealized model molecular cloud. In a subsequent paper we consider the effects of the grain fluid(s).

Langer (1978) studied the stability of a model molecular cloud (infinite in extent and uniform in density and magnetic field) with respect to small-amplitude, adiabatic perturbations in the presence of ambipolar diffusion. For propagation along the magnetic field lines, he recovered, as one would expect intuitively, the Jeans dispersion relation and instability in the absence of the magnetic field. He then investigated the wave propagation perpendicular to the field lines. He showed that the Jeans instability is still present, that the critical wavenumber for instability is independent of the magnetic field strength, but that the growth rate depends on both the field strength and the degree of ionisation. Aside from two spurious curves in his Figure 1, which exhibits the growth rate and decay time of some modes, our results for propagation of the low-frequency modes perpendicular to the magnetic field are in agreement with Langer’s – he ignored the high-frequency ion modes. Yet even in this, previously studied case, we offer new analytical expressions and new physical insight and interpretation of the results. Moreover, we present not only the eigenvalues (frequencies or, equivalently, phase velocities) as functions of wavelength but the eigenvectors as well (i.e., material velocities, densities, magnetic-field components, etc.). We also study propagation at arbitrary angles with respect to the magnetic field, and we offer a thorough discussion of the wave modes, not just the ambipolar-diffusion–induced instability.

Pudritz (1990) revisited Langer’s problem (with the minor difference of considering isothermal perturbations) but introduced a new effect: he assumed that there exists a power-law spectrum of small-amplitude waves, and then he studied the effect that this spectrum has on the ambipolar-diffusion–induced, Jeans-like instability. 2 He concluded that the slope of the spectrum (considered as a function of wavelength) has an important effect on the growth rate of the instability; the steeper the spectrum, the greater the growth rate. (The growth rate of gravitationally-driven ambipolar diffusion, however, cannot possibly exceed the free-fall rate.)

Several other papers have appeared in print since 1990, studying different aspects of weakly ionised systems, focusing usually on the stability of certain MHD modes or shocks, especially as it may relate to the formation of structures in molecular clouds and/or on the effect of the grain fluid(s) on the allowable wave modes or shocks (e.g., Wardle 1990; Balsara 1996; Zweibel 1998; Kamaya & Nishi 1998, 2000; Mamun & Shukla 2001; Cramer et al. 2001; Falle & Hartquist 2002; Tytarenko et al. 2002; Zweibel 2002; Ciolek et al. 2004; Lim et al. 2005; Oishi & Mac Low 2006; Roberge & Ciolek 2007; van Loo et al. 2008; Li & Houde 2008). In this paper we present a general theory of the propagation, dissipation and growth of MHD waves in partially ionised media in three dimensions, with emphasis on mathematical transparency of the formulation and analytical solution of the problem, the physical understanding and interpretation of all modes, including their eigenvectors, the many critical wavelengths that exist and which separate regimes dominated by different waves or instabilities, and on specific features relevant to the evolution of molecular clouds. As mentioned above, even when a particular result agrees with previous work, we offer new insight into its physical understanding.

In § 2 we present the equations governing the behaviour of a weakly ionised, magnetic, self-gravitating interstellar cloud. The equations are linearised, Fourier-analyzed, and put in dimensionless form. The free parameters of the problem are identified, their physical meaning explained, and their typical values given. The different hydromagnetic modes and their dependence on wavelength for different directions of propagation relative to the magnetic field are calculated and explained physically in § 3. Analytical expressions for the phase velocities, damping timescales, growth timescales, including critical or cutoff wavelengths, are also obtained. A physical discussion of the eigenvectors is an integral part of this presentation. Section 4 summarizes some of the results and their relevance to the formation of protostellar fragments (or cores) and to other observable phenomena. It also gives in two Tables all the critical wavelengths and the ranges of wavelengths in which different modes can exist in molecular clouds, for propagation parallel, perpendicular, and at arbitrary angles with respect to the magnetic field.

2 Formulation of the Problem

2.1 Basic Equations

We consider a weakly ionised medium (e.g., an interstellar molecular cloud) consisting of neutral particles ( with a helium abundance by number; subscript n), electrons, and singly-charged positive ions (subscript i). For specificity we assume that the ions are molecular ions (such as ); for the densities of interest in this paper (), this is sufficient since atomic ions (such as or ) are less abundant (a result of depletion of metals in dense clouds) and, in any case, they have masses comparable to that of (for more detailed treatments of the chemistry, see Ciolek & Mouschovias 1995, 1998, or the appendix of Mouschovias & Ciolek 1999). Interstellar grains, which have been shown to have significant effects on the formation and contraction of protostellar cores (Ciolek & Mouschovias 1993, 1994, 1995) and in the opaque phase of star formation (Tassis & Mouschovias 2007a, b, c; Kunz & Mouschovias 2009, 2010) are neglected in this analysis; they are accounted for in a subsequent paper.

The magnetohydrodynamic (MHD) equations governing the evolution of the above two-fluid system are


where and are, respectively, the density and velocity of species ,


is the time-derivative comoving with species , the gravitational potential, the neutral pressure, the magnetic field, the temperature, and the heating and cooling rates (per unit volume) of the neutral gas, and


the neutral-ion and ion-neutral mean collision (i.e., momentum-exchange) times. The quantity is the universal gravitational constant, and is Boltzmann’s constant; is the mean mass of a neutral particle in units of the atomic-hydrogen mass and is equal to 2.33 for a gas with a helium abundance by number. The quantities and in the ion mass continuity equation (1b) are, respectively, the cosmic-ray ionisation rate and the coefficient for dissociative recombination of molecular ions and electrons (in ). In writing equation (1b), we have used the condition of local charge neutrality , where and are, respectively, the number densities of ions and electrons. We may use the assumption of local charge neutrality because the various HM modes of interest here have frequencies much smaller than the electron plasma frequency (where is the abundance of electrons relative to the neutrals, and is equal to the degree of ionisation for weakly-ionised systems); hence, any excess charge density is quickly shielded by the mobile electrons, so that for timescales . Since we neglect the effects of grains in this paper, capture of ions onto grains is not included on the right-hand side of equation (1b) as a sink term for ions.

Heat conduction and viscosity are not important for the densities and lengthscales of interest (, ), and are therefore ignored as possible sources of heating/cooling in the model clouds. (Note that the left-hand side of equation [1e] is equal to , where is the entropy per gram of matter.)

Because , we include only the neutral density as a source term in Poisson’s equation (1h). Similarly, the gravitational and thermal-pressure forces (per unit volume) on the plasma (ions and electrons) have been neglected in the plasma force equation (1d). One can easily show that, for the physical conditions in typical molecular clouds, they are completely negligible in comparison to the magnetic force exerted on the plasma, except in a direction almost exactly parallel to the magnetic field. Ignoring them parallel to the magnetic field implies that we are neglecting the ion acoustic waves and the (extremely long-wavelength) Jeans instability in the ions.

The quantity in equations (3a) and (3b) is the elastic collision rate between ions and neutrals. For collisions, (McDaniel & Mason 1973). The factor 1.4 in equations (3a) and (3b) accounts for the inertial effect of He on the motion of the neutrals (for a discussion, see § 2.1 of Mouschovias & Ciolek 1999).

2.2 Linear System

To investigate the propagation, dissipation, and growth of HM waves in molecular clouds, we follow the original analysis by Jeans (1928; see also Spitzer 1978, § 13.3a; and Binney & Tremaine 1987, § 5.1), and assume that the zeroth-order state is uniform, static (i.e., ), and in equilibrium. 3 We consider only adiabatic perturbations; therefore, the net heating rate () on the right-hand side of equation (1e) vanishes.

We write any scalar quantity or component of a vector in the form , where refers to the zeroth-order state, and the first-order quantity satisfies the condition . We thus obtain from equations (1a) - (1i) the linearised system


Equation (4b) has been simplified by using the relation


which expresses equilibrium of the ion density in the zeroth-order state, as a result of balance between the rate of creation of ions from ionisation of neutral matter by high-energy () cosmic rays and the rate of destruction of ions by electronmolecular-ion dissociative recombinations. This relation allowed us to replace by in equation (4b). The quantity in equation (4b) is the degree of ionisation (where and are the number densities of ions and neutrals in the unperturbed state). For an ideal gas (with only translational degrees of freedom), in equation (4e).

We seek plane-wave solutions of the form , where is the propagation vector, the frequency, and the amplitude (in general, complex) of the perturbation. Equations (4a) - (4i) reduce to




is the adiabatic speed of sound in the neutrals, and


is the (one-dimensional) neutral free-fall timescale. The quantities , , and have been eliminated by using equations (4e), (4g), and (4h), respectively.

2.3 The Dimensionless Problem

We put equations (6a) - (6f) in dimensionless form by adopting , , , and as units of density, magnetic field strength, time, and speed, respectively. The implied unit of length is , which is proportional to the one-dimensional thermal Jeans lengthscale (see § 3.1.1). For convenience, we adopt a cartesian coordinate system such that the propagation vector is in the -direction and the zeroth-order magnetic field is in the -plane at an angle with respect to (see Fig. 1). Then the three unit vectors are

Figure 1: Coordinate system used in analysing the linearised hydromagnetic equations. The -axis is aligned with the propagation vector , and the magnetic field is in the -plane at an angle with respect to (see eqs. [9a] - [9c]).

One may write , and the dimensionless form of equations (6a) - (6f) can be written in component form as


Note that equation (10l) is redundant in that it gives the same information as equation (10e), namely, that there cannot be a nonvanishing component of the perturbed magnetic field in the direction of propagation.

The dimensionless free parameters appearing in equations (10a) - (10l) are given by

they represent, respectively, the neutral-ion collision time, the ion-neutral collision time, the ion Alfvén speed , and the electronmolecular-ion dissociative recombination rate per unit ion mass. In evaluating the numerical constants in equations (11a) - (11d) we have used and ; we have also normalized the ion mass to that of (= 29 amu). For any given ion mass and mean mass per neutral particle (in units of ), the ion mass fraction and the cosmic-ray ionisation rate are not free parameters in the problem; the former is determined by the ratio and the latter by the product , where is the dimensionless dissociative-recombination coefficient for molecular ions.

We note that , where is the collapse retardation factor, which is a parameter that measures the effectiveness with which magnetic forces are transmitted to the neutrals via neutral-ion collisions (Mouschovias 1982), and appears naturally in the timescale for the formation of protostellar cores by ambipolar diffusion (e.g., see reviews by Mouschovias 1987a, § 2.2.5; 1987b, § 3.4; 1991b, § 2.3.1; and discussions in Fiedler & Mouschovias 1992, 1993; Ciolek & Mouschovias 1993, 1994, 1995; Basu & Mouschovias 1994, 1995). It is essentially the factor by which ambipolar diffusion in a magnetically supported cloud retards the formation and contraction of a protostellar fragment (or core) relative to free fall up to the stage at which the mass-to-flux ratio exceeds the critical value for collapse. It is discussed further in § 3.2.1.

Equations (10a) - (10k) govern the behaviour of small-amplitude disturbances in a weakly ionised cloud; they (without eq. [10e]) form a homogeneous system. In general, the dispersion relation can be obtained by setting the determinant of the coefficients equal to zero. To each root (eigenvalue) of the dispersion relation there corresponds an eigenvector (or “mode”), whose components are the dependent variables appearing in equations (10a) - (10k). (Note that, once the dependent variables in eqs. [10a] - [10k] are known, one may use eqs. [4e], [4g], and [4h] to solve for the perturbed quantities , , and , respectively.) Since, in general, is complex, modes with decay and those with grow exponentially in time. In what follows we investigate the propagation, dissipation, and growth of the allowable HM modes in typical interstellar molecular clouds.

3 Solution, Physical Interpretation, and Applications

For specificity, we consider a representative molecular cloud of density , magnetic field strength , temperature , dissociative recombination rate , and cosmic-ray ionisation rate , implying a degree of ionisation (and, hence, ion mass fraction ). The unit of time (see eq. [8]) for this model is equal to , and the unit of speed is (see eq. [7]). Hence, the unit of length is . The four (dimensionless) free parameters of the problem (see eqs. [11a] - [11d]) are: the ion-neutral collision time , the neutral-ion collision time , the Alfvén speed in the ions , and the dissociative-recombination coefficient . (These imply a dimensionless cosmic-ray ionisation rate .)

In the following subsections we present the solutions for propagation along (), perpendicular (), and at intermediate angles (, , and ) with respect to the unperturbed magnetic field (see Fig. 1). We use the velocity vector in relation to as defining the polarisation of each wave mode. Modes that have only , are said to be longitudinally polarised; it follows from equation (10a) that these modes are compressible, i.e., . Modes that have , are said to be transversely polarised; they are incompressible, i.e., . Note in what follows that, since the thermal pressure in the plasma has been neglected, no ion sound waves are present.

3.1 Propagation Along ()

For , equations [10a] - [10k] become uncoupled in the three mutually orthogonal directions , , and . The modes polarised in the -direction are given by (see eqs. [10a] - [10d])


Equations (10f) - (10h) yield for the modes with motions in the -direction


while equations (10i) - (10k) for the modes polarised in the -direction become


Longitudinal Modes for : Eigenfrequencies and Eigenvectors

From equations (12a) - (12d) the dispersion relation for the longitudinal modes is


The four eigenvalues and eigenvectors as functions of dimensionless wavelength () are obtained from direct numerical solution of the eigensystem (12a) -(12d) and displayed in Figures 2 and 3. Figure shows the absolute value of the phase velocity (, where ). The damping (or dissipation) timescale (, , where ) is exhibited in Figure ; growth times (, ) are shown in Figure . The absolute values of the -components of the eigenvectors for the various modes are shown in Figure , while those of , , and are displayed in Figures , , and , respectively. Note that, in Figures - , all the eigenvectors have been normalised to unity; i.e., they satisfy the condition .

Figure 2: Eigenvalues of longitudinal modes as functions of wavelength, normalised to , at an angle of propagation with respect to . Because of degeneracy, there are a total of four different modes at each wavelength. Each curve is identified by a label (a mnemonic for the mode it represents) as explained in Table 1. (a) Absolute value of the phase velocity , normalised to . Also overplotted as boxes with interior ’s are values obtained from eq. (19). (b) Damping timescale , in units of . Also shown (open circles) are values calculated from eq. (21). (c) Growth timescale , normalised to . Also plotted (open circles) are values calculated using eq. (21).
Figure 3: Magnitudes of eigenvectors of longitudinal modes for as functions of wavelength, normalised as in Fig. 2. As in Fig. 2, there are a total of four modes. All eigenvectors have been normalised to unity. Each curve is identified by a label (a mnemonic for the mode it represents) as explained in Table 1. (a) Longitudinal (-) component of the neutral velocity, . (b) Longitudinal (-) component of the ion velocity, . (c) Neutral density . (d) Ion density .

The first mode for this system of equations is a nonpropagating, ion collisional-decay mode, i.e., the ions are streaming through a sea of fixed neutrals. Their motion decays because of ion-neutral collisions. It is characterized by . Solving equations (12a) - (12d) under these conditions (or, equivalently, solving eq. [15] in the limit ), one finds


Hence, and (see Fig. 2, line labeled “i,coll”). This mode is independent of wavelength (see Figs. , , lines labeled “i,coll”).

The second mode is one in which density enhancements in the ions rapidly decay by dissociative recombinations of molecular ions and electrons; because the degree of ionisation is so small (), this mode does not involve any motion of the neutrals and leaves the neutral density essentially unchanged (see Figs. and , lines labeled “i,rec”). Solving equation (12b) with and , one finds that and


which is equal to for the model described here. It is again the case that this mode is independent of (see Figs. and - , lines labeled “i,rec”).

The remaining two longitudinal modes are low-frequency modes, with . Because the inertia of the ions along the magnetic field is small (), the neutrals are able to sweep up the ions, and, as a result, (see Figs. 3 and 3). For these conditions, equation (15) yields the thermal Jeans modes


(e.g., see Chandrasekhar 1961, Ch. XIII; or Spitzer 1978, § 13.3a). Therefore, for ,


i.e., the two acoustic waves have the same phase velocity, modified by gravity, but propagate in opposite directions (along the field lines). In the limit , and ; i.e., these modes are undamped sound waves (recall that the unit of speed is ). At longer wavelengths, gravitational forces become increasingly more important and the phase velocity becomes less than unity (see Fig. ). The waves are gravitationally suppressed (i.e., ) at wavelengths greater than the thermal Jeans wavelength


(, dimensionally). For (i.e, ), it follows from equation (18) that each of the Jeans modes splits into two separate, conjugate modes. One is a gravitational growth (or fragmentation) mode, with timescale


(see Fig. ). This is the classical Jeans instability. As , ; dimensionally, this is just the free-fall timescale, . The corresponding eigenvector is labeled as “ff+” in Figures - . The other mode is one of exponential decay, with damping timescale also given by equation (21) (see eq. [18]); it is the curve labeled by “ff” in Figure . The eigenvector, also labeled by “ff”, is shown in Figures - . This mode is one in which an initial density enhancement causes expansive motion, opposed by gravity, at such a rate that the density enhancement decreases to zero at the same time that the velocity vanishes. Hence, this is a monotonically decaying mode; no wave motion is involved. It is similar to the well-known classical cosmological problem of an expanding “flat” universe. Note that, as , .

In order to better understand the neutral thermal (Jeans) modes, we examine more closely the eigenvectors (and their features shown in Figures - ). We substitute equation (18) in equation (12a), and we use the normalisation condition