### i s s n 0 0 6 5 - 3 7 1 3

**I N S T I T U T O ' A E R O N O M I E S P A T I A L E O E B E L G I Q U E ** 3 - A v e n u e C i r c u l a i r e

### B • 1180 B R U X E L L E S

**A E R O N O M I C A A C T A **

### A - N° 354 - 1990

### Nonlinear supersonic motions of magnetized plasma flow with mass loading

### by

### I. KH. KHABIBRAKHMANOV and F. VERHEEST

**B E L G I S C H I N S T I T U U T V O O R R U I M T E - A E R O N O M I E **

3 - flmglaan
**F O R E W O R D **

### This paper was m a d e possible by the cooperation between the BISA and the Institute of Space Research (Moscow, USSR) and will be published in the Journal of *Geophysical Research (Space Physics) *

**V O O R W O O R D **

### Dit artikel is mogelijk gemaakt door de samenwerking tussen het BIRA en het Instituut voor Ruimte-onderzoek (Moskou, USSR) en zal in het Journal of Geo- *physical Research (Space Physics) gepubliceerd worden *

**A V A N T - P R O P O S **

### Cet article a été rendu possible grâce à la coopération entre l'IASB et l'Institut de Recherches Spatiales (Moscou, URSS) et sera publié dans le Journal of Geophysical *Research (Space Physics) *

**V O R W O R T **

### Diese Veröffentlichung ist dank der Zusammenarbeit zwischen dem BIRA

### und dem Raumforschungsinstitut (Moskou, USSR) ermöglicht worden u n d wird in

### Nonlinear supersonic motions of magnetized plasma flow with mass loading

### by

**I. Kh. Khabibrakhmanov and F. Verheest **

**Space Research Institute, USSR Academy of Sciences * **
**Profsoyuznaya 8 4 / 3 2 , Moscow 117810, USSR **

**Instituut voor theoretische mechanika, Rijksuniversiteit G e n t * * **
**Krijgslaan 281, B - 9 0 0 0 Gent, Belgium **

**A b s t r a c t **

### A kinetic description of magnetized plasma flow with mass loading is devel- oped and leads to a nonlinear wave equation for small scale dispersive magnetosonic plasma motions superimposed on large scale motions due to mass loading. The mass loading can drive quasistationary, finite-amplitude magnetosonic type waves in a su- personic plasma flow, even though small amplitude linear disturbances are damped, but only when the magnetosonic Mach number drops to 2 during the slowing down of the mass ioaded plasma flow. Nonlinear motions with constant amplitude could be permanent features of a mass loaded plasma flow, as in the neighborhood of comets and planets \yith extended atmospheres but negligible magnetic fields. The relation between the wave amplitude and the local magnetosonic Mach number is computed and can easily be compared with in situ measurements in cometary plasma environ- ments.

**S a m e n v a t t i n g **

### Een kinetische beschrijving is ontwikkeld voor een gemagnetizeerd stromend

### plasma met massalading en leidt tot een niet-lineaire golfvergelijking voor klein-

### schalige dispersieve magnetosonische plasmabewegingen, gesuperponeerd op groot-

### schalige bewegingen te wijten aan de massalading. Deze massalading kan quasista-

### tionaire magnetosonische golven met eindige amplitude opwekken in een supersonisch

### stromend plasma, zelfs al zijn lineaire golven met kleine amplitude gedempt, maar

### alleen als het magnetosonisch Machgetal naar 2 daalt gedurende de vertraging van

### het stromend plasma door de massalading. Niet-lineaire golven met konstante am-

### plitude zouden dus blijvende kenmerken kunnen zijn van stromende plasma's met

### massalading, zoals in de omgeving van kometen en planeten met een uitgestrekte

### atmosfeer maar met verwaarloosbaar magneetveld. Het verband tussen de golfam-

### plitude and het lokaal magnetosonisch Machgetal is berekend en kan gemakkelijk

### vergeleken worden met in situ metingen in kometaire plasma's.

**R é s u m é **

### On a développé une description cinétique du flux d'un plasma magnétisé avec chargement de masse, résultant en une équation nonlinéaire pour des mouve- ments dispersifs et magnétosonores à petite échelle, superposés à des mouvements du plasma à grande échelle dûs au chargement de masse. Ce chargement de masse peut générer des ondes magnétosonores quasistationnaires d'amplitude finie dans un flux de plasma supersonique, bien que des perturbations de petite amplitude soient amor- ties, mais seulement quand le nombre de Mach magnétosonore tombe à 2 pendant la décélération du flux de plasma avec chargement de masse. Ces ondes nonlinéaires à amplitude constante pourraient être des caractéristiques permanentes d'un flux de plasma avec chargement de masse, comme au voisinage des comètes et planètes avec une atmosphère étendue mais un champ magnétique négligeable. La relation entre l'amplitude de l'onde et le nombre de Mach magnétosonore local est calculée et peut facilement être comparée aux mesures faites in situ dans l'environnement des plasmas cométaires.

**Z u s a m m e n f a s s u n g **

### Eine kinetische Beschreibung eines magnetisierten Plasmaflusses mit Massen-

### ladung ist entwickelt und führt zu einer nichtlinearen Wellengleichung für disper-

### sive, magnetosonore Plasmabewegungen auf kleinen Maßstab, superponiert auf Be-

### wegungen mit großem Maßstab wegen der Massenladung. Diese Massenladung kann

### quasistationären, magnetosonoren Wellen mit konstanter Amplitude antreiben in

### einem Überschallfluß , eben wenn lineare, infinitesimale Störungen gedämpft wer-

### den, aber nur wenn die magnetosonore Machzahl zu 2 senkt während der Plasmafluß

### mit Massenladung langsamer wird. Solche nichtlineare Bewegungen mit konstanter

### Amplitude könnten permanente Kennzeichen des Plasmaflusses mit Massenladung

### darstellen, wie in der Nähe von Kometen und Planeten mit ausgedehnter Atmo-

### sphäre aber vernachlässigbarem Magnetfeld. Die Relation zwischen die Wellenampli-

### tude und die lokale magnetosonore Machzahl ist berechnet worden und kann leicht

### mit in situ Meßwerten in kometaren Plasmen verglichen worden.

**1. I N T R O D U C T I O N **

### There are various astrophysical situations where a magnetized plasma flows through a neutral gas. One prominent example is the interaction of the solar wind with extended cometary atmospheres, studied in detail following the recent success- ful cometary missions (see review by Galeev, 1987), or with planets having relatively extended atmospheres but small or no intrinsic magnetic fields. Among the more distant objects one could mention stellar winds, or astrophysical plasma jets from the nuclei of active galaxies, or again expanding shells of supernovae remnants pass- ing through and interacting with the neutral component of the interstellar medium *[Petelski et al., 1980; Petelski, 1981]. Inside the Jovian magnetosphere there is yet * another example of such an interaction: the Io plasma torus resulting from the ion- ization of material evaporating from the surface of Io [Galeev and Khabtbrakhmanov, 1986; Horton and Smith, i988j.

### The interaction of a plasma flow with a neutral gas gives rise to a so-called mass loading. By this one means that newly ionized particles, created by several ionization processes such as photoionization, electron impact or dissociative ioniza- tion, are incorporated into the main plasma flow. The relative motions between the original plasma and the newly ionized particles induce electric fields, which accelerate the new ions and give them drift .velocities Vp = E x B / B

*2*

*. There is, of course, the * feedback effect on the original plasma flow, which is decelerated so as to conserve momentum and energy in the combined system. In the case of slow electromagnetic field variations (adiabatic motions of the particles), the new plasma particles acquire a thermal velocity equal to the drift velocity, due to their motion transverse to the magnetic field.

### So the main effect of the mass loading on the plasma flow is its deceleration,

### although the accompanying pressure gradients can in some cases compensate the

### inertial forces and ensure a steady state situation. This was e.g. in the case of the

### solar wind interaction with comets originally considered by Walks [1973] and Wallis

*and Ong [1975] and more recently by Galeev, Cravens and Gombosi [1985]. They *

### found that smooth deceleration of the supersonic solar wind flow is possible until the

### mass flux of the flow reaches a critical value (pu)

^{c}### = 4/3, corresponding to a local

### Mach number M = 1. This means that a stationary picture cannot be the whole

### story and t h a t a shock transition has to occur somewhere before the sonic point of the flow.

### In the case of planets with intrinsic magnetic fields the interaction region between the solar wind and the magnetosphere is usually very small compared to the size of planet. This means t h a t the braking of the solar wind flow occurs in a short distance and the characteristics of the stationary bow shock are determined almost exclusively by the properties of the solar wind far u p s t r e a m . In the case of comets, however, having no magnetic field of any significance, no rigid obstacle is present to play the role of a piston generating a shock. Physically this means t h a t the interaction region is very wide and its size is determined by the coma. Nevertheless the solar wind flow possesses another intrinsic spatial scale, the dispersive length scale, which was the only i m p o r t a n t one in the case of a magnetized planet. It is precisely the interaction of these two very different scales which determines the bow shock in the case of comets.

### T h e ratio of the dispersive scale to the characteristic size of the coma intro- duces in a n a t u r a l way a small p a r a m e t e r and one can then exploit the weak inter- action between plasma motions on such vastly different scales. T h a t is, small scale motions, on scales much shorter t h a n the ionization scalelengths, have to influence the plasma flow. In a hot collisionless plasma the low frequency magnetosonic type wave motions exert a m a j o r influence, higher frequency disturbances being globally of much less of importance.

*« Galeev and Khabibrakhmanov [1989] considered the collisonless interaction * of small scale dispersive motions with large scale motions induced by mass loading, in the case of solar wind interaction with cometary atmospheres. This interaction leads to finite amplitude magnetosonic modes in the region of mass loaded plasma flow with magnetosonic Mach numbers less t h a n 2. This Mach number is defined by characteristics of the large scale motions. It is thus n a t u r a l to suppose t h a t the real cometary shock must be somewhere in this region, t h a t is, the interaction of large and small scale motions is indeed responsible for the bowshock formation.

### Numerical simulations [Galeev, Lipatov and Sagdeev, 1985; Omidi and Win-

*ske, 1986] show t h a t the cometary bowshock Mach number is close to 1.7 ~ 2 and *

### is only weakly dependent on the supply rate of neutral material from the cometary

### nucleus. Another consequence of this interaction between motions on different scales is the damping or growth of nonlinear waves. As a result, there are quasistationary nonlinear waves, propagating with almost constant amplitudes.

### In the present paper we will analyze quas is tat ion ary nonlinear motions in general, following the approach outlined by Galeev and Khabibrakhmanov [1989].

### Both that paper and the present one run along similar lines, but here we fully take all nonlinear and mass loading effects into account, contrary to previous treatments. We expand the solution of the Vlasov equation for the plasma flow particles (Appendix C) in appropriate variables // and <f>, where (J. = v\/2B is the magnetic moment and 4> the modified gyrophase angle of a particle (Appendix B). The right hand side of this Vlasov equation is determined by the ionization rate of the neutral particles (Appendix A). The expansion of the solution allows us to calculate the total current.

### Upon substitution of the resulting currents for each plasma species into Maxwell's equations we get a set of equations describing the motion of the magnetized plasma flow. The analytical analysis is then straightforward, as the equations are already conveniently expanded in the derivatives of the electromagnetic fields. So in Section 2 we present the kinetic description of slow motions of the mass loaded plasma flow. In Section 3 we consider the ionization scale motions of the flow and in Section 4 motions on a much shorter dispersive scale. Such a division is convenient and justified here because of the very large differences between these scales. Finally in Section 5, the interaction between motions on these two different scales determines the evolution of the quasistationary nonlinear waves.

**2. K I N E T I C D E S C R I P T I O N OF M A G N E T I Z E D , M A S S L O A D E D ** **P L A S M A F L O W **

### For simplicity we will consider the motion of a mass loaded plasma flow only

### in the plane transverse to the magnetic field direction (2-axis), magnetic field which

### has thus only one component. The Vlasov equation for the distribution function of

### the plasma particles in the variables n = v\/2B (the magnetic moment) and 4> (the

### modified gyrophase angle, see Appendix A) is given in Appendix B. The newly ionized

### particles, embedded now in the main plasma flow, are considered to be initially cold,

### as the thermal speed of the neutral particles is very small compared to the local

### drift velocity of the plasma. This assumption yields a relatively simple form of the

### distribution function of newly ionized particles (A.l) and of the source term in the Vlasov equation for the plasma species (A.4). We consider quasistationary motions of the plasma with spatial scales much greater , than the ion gyroradius r

^{c}

### , and the electron gyroradius r

^{c e}

### . Using this approximation we expand the solution of the Vlasov equation for the ions in the small parameter e = 0(w/f2) = 0{kr

^{c}*i), where 0 * is the gyrofrequency of the ions and u and k are the inverse time and spatial scales of the plasma motion, respectively. The expression for the solution of the Vlasov equation up to the third order in e is given in Appendix C. With the help of this solution one calculates the total current j :

### j

^{ =}

### / v F

^{Q}

### ( r , v , 0 d

^{3}

### v , (2.1)

### a

### where the subscript a denotes the plasma species and e

*a*

### is the charge of particles of species a . This expression can be used to obtain the equations of motion for the plasma flow from Maxwell's equation

### c

2### V x B = ^ + - j . (2.2)

*Ot £ o *

### The y-component of this is the equation of motion for the plasma along the x-axis.

### The drift velocity along the x-axis, vox = ƒ = E

^{y}*/B is in a first approximation the * bulk velocity of the plasma.

### Now we will consider only one-dimensional motions along the x-axis. This implies that the drift velocity v

^{Dy}* = -g = -E*

^{x}*/B along the y-axis is zero, although * it has nonvanishing gradients. It is then easy to determine the evolution of g from the x-component of equation (2.2):

*3P ~ n*

^{ r }### 2 dx mB

^{2}### 2 B

*\fTB ~ v *

### + - T— = 0 . (2.3)

### Here we used up to second-order terms (the orders of the terms refer to powers of

### l/f2), because g occurs only in the second or higher order terms in the equation for

### the motion along the stagnation line:

**dEr, o dB eB**

### f P 2P _ j//

### 3 P

**+**** n U****2**** 2mB****2 **

### 1 / I I P

U### T

3### f

3 1 1~ Ô^{2 B x +}**B**^{ f} **~ **

**B**

**B**

*2*

**f d f**** P** a

**P**

^{3}

### n + r ^ ^ + B

### I I I

### 2ft

2** I dx mB**

2**I dx mB**

### dx

3### m

2### £

3### y

### _ f f

2### - , / ! ! - — — - f ^ - )

3### '

2### ' ) } = o .

**U**

^{2}

**\ B\ 2B 3 dx**

^{2}

**B\2Bj J j**### (2.4) **In equations (2.3-4) the operator T is denned as T = d/dt + fdjdx — gd/dy. For ** **the sake of simplicity, all variables will eventually be supposed independent of y, so ** that later we only consider one-dimensional motions of the plasma flow. Also, the **standard designation for derivatives d./dx = .**

**In equations (2.3-4) the operator T is denned as T = d/dt + fdjdx — gd/dy. For**

**the sake of simplicity, all variables will eventually be supposed independent of y, so**

**standard designation for derivatives d./dx = .**

**x**### is used.

**We have to add equations governing the evolution of the plasma density n, ** r oo

**We have to add equations governing the evolution of the plasma density n,**

**n = B Fo(r,M)dM, (2.5) **

**n = B Fo(r,M)dM, (2.5)**

**Jo**### of the ion pressure P ,

### P f°°

**- = B**

**- = B**

^{2}### F

^{0}

### (r,Ax)/xdAi, (2.6)

**m Jo **

### of the second moment II (in the magnetic moment y, of the particles),

TT f ° °

### = B

3### P

0### (r,M)M

2### dM, (2-7)

**m****i**** J****0 **

**and of the magnetic field B. However, it is easier to use equations governing the **

**and of the magnetic field B. However, it is easier to use equations governing the**

### changes in n / B , P / B

^{2}

### and I I / B

^{3}

### , because of their invariance in the case of zero

### production rate u. This is a consequence of the adiabatic conservation of the magnetic

### moment of the particles during slow motions:

**For the evolution of the magnetic field we can use Maxwell's equation V x E + d B / d t = ** 0, which in our notations is written as

*TB = (g*

^{y}*-f*

^{x}*)B. (2.11) *

### So we have a complete set of equations (2.3-4),(2.8-11) describing the slow motions of the plasma, As can be seen from the method of expansion of the distri- bution function in Appendix C, slow motions here mean motions with characteristic time and space scales much greater than the ion gyroperiod and gyroradius, respec- tively.

**3. QUASISTATIONARY PLASMA MOTIONS: **

**LARGE SCALE MOTIONS **

### Let us now consider the quasistationary motions described by the set (2.3- 4). As was stated in the Introduction, we take into account motions simultaneously on two very different spatial scales and it is hence advantageous to distinguish con- tributions on different scales from the very outset. Every variable of the system is supposed to have two parts, one purely stationary and another one with a weak time dependence. The former describes the large space scale variations of the plasma flow, whereas the latter takes dispersive effects into account. For example, the flow veloc- *ity ƒ along the stagnation line has a stationary part f(x) and a part f(t,x) which * depends only weakly on time.

### Neglecting the small scale motions in the set (2.3-4),(2.8-11) yields the

### following magnetohydrodynamic equations:

*ƒ d h v * *dx B B' * *d_P_*

*=*

* vrn] *

*dx B*

*2*

* 2B*

*2*

* ' * *d_U_ _ um*

*2*

*f*

*3*

*dx B*

*3*

* ~ 4B*

*z*

* ' * *dBf _ *

*dx *

### T h e large scale motion is precisely the one used by Galeev, Craven and Gombosi [1985] to describe the solar wind flow in the vicinity of comets. We can use their solutions for (3.1), in the limit of small magnetic field pressure (neglecting it in the first equation):

*P = mnoo/^ - mhf*

*2*

*, *

~ r r I ^ / o o

* foo *

*n j — ttoo / o o *

### 3ƒ 3 / 2 (3.2)

^ o o / o o ^ o o / o o 4 7 r V gr n00/ o o r

### Here the subscript oo denotes the undisturbed solar wind values. T h e plasma pro- duction rate u depends on the production rate Q

*n*

### of neutral cometary gas, on a constant gas outflow velocity V

*g*

*, on a characteristic time T for photoionization a n d * on the radial distance r f r o m the cometary nucleus, hence u = Q

^{n}

### / ( 4 7 r V ^ r r

^{2}

### ) .

### It can now be seen t h a t the characteristic space scale for large scale plasma motions is

**RL**

### = .

u 9 n### , • (3-3)

### T h e expressions (3.2) completely define these large scale motions, except near the

### point with large gradients where the local Maçh number, defined through M

2### =

*mhf*

*2*

*/2P = 2(ƒ - 1 / 4 ) / ( 1 - ƒ), goes to 1. This means ƒ = 1/2, with ƒ =*

/ ' / / o o -
**4. QUASISTATIONARY PLASMA MOTIONS: **

**DISPERSIVE SCALE MOTIONS **

### Now we will take into account the dispersive effects on the mass loaded solar wind flow. This is easily done by taking advantage of the fact that the ratio of the dispersive scale ao to the plasma production scale Re (defined in (3.3)) is so small. The spatial derivatives of the large scale motions will then be smaller than those of the small scale motions by the same ratio. The following analysis is in fact based on the expansion of the set (2.3-4),(2.8-11) in the small parameter

**CID/RL-**### Before doing that, however, we note that in the opposite limit (Re = 0, signifying the absence of mass loading) the set (2.3-4),(2.8-11) describes nonlinear magnetosonic **type motions (see e.g. Mikhailovskii and Smolyakov [1985]). **

**type motions (see e.g. Mikhailovskii and Smolyakov [1985]).**

**For quasistationary motions, where in the operator T the time derivatives ** **d/dt are much smaller than the fd/dx terms, one reduces the set (2.3-4),(2.8-11) to ** one equation for perturbations in the magnetic field. These perturbations, denoted **by b and defined through B/B = 1 + b(t,x — ut), obey a KdV-equation **

**For quasistationary motions, where in the operator T the time derivatives**

**d/dt are much smaller than the fd/dx terms, one reduces the set (2.3-4),(2.8-11) to**

**by b and defined through B/B = 1 + b(t,x — ut), obey a KdV-equation**

### 2P + £

2### /mo \

b### ^

+### ^

+ a### 2

D b x x x =### o. (4.1)

*P(f-u) J *

### Here p = mh is the mass density of the plasma and ap the dispersion length for the magnetosonic wave, defined through

**a**

**a**

**D**** = ** (4.2)

**=**

### We see that in the absence of mass loading (u = 0, that is, the values of ri/B, P / B

^{2 }

**and Tl/B**

**and Tl/B**

**3**### are adiabatically conserved) the set (2.3-4),(2.8-11) describes soliton-type

### nonlinear motions: 6(0)

### 6 = 7 Y T 1 777' (4-3) **cosh («(x — ut J) **

**cosh («(x — ut J)**

### where

### ( ƒ - " )

### 1 . I V~

2*V*

=
### are the soliton amplitude, its characteristic width and the velocity of the linear mag- netosonic waves, respectively.

### T h e preceeding results are well known. They include large scale motions, describing the smooth deceleration of the solar wind flow, as well as soliton-type wave motions on the small scale (dispersive motions). It is worth noting t h a t the description used in this paper is purely a kinetic one and thus automatically obviates problems concerning the precise calculation of the viscosity tensor in a h y d r o d y n a m i c t r e a t m e n t (see Kennel and Sagdeev [1967]; Macmahon [1968]; Kennel [1968]). T h e kinetic description enables us to reproduce results of Mikhailovskii and Smolyakov [1985] derived by another approach — they used linear kinetic theory to calculate thé dispersive terms and then a hydrodynamic model to describe nonlinear effects in the wave equation. Moreover, in principle only a full kinetic t r e a t m e n t can give an accurate c o m p u t a t i o n of the effects of different w a r m plasma components on the nonlinear motions, in particular on the dispersive lengthscale. Although Appendix C has all the expressions needed for such calculations, it is not the aim of the present paper to pursue this aspect f u r t h e r .

### It is also well known t h a t if we include small dissipative terms, the soliton solution has two different kinds of evolution, depending on how the dissipation is cat- alogued. If the dissipative t e r m is proportional to b ("hydrodynamic" dissipation), we have the motion of a soliton with slowly varying amplitude. In this case there is no shock. If on the other h a n d the dissipative t e r m is proportional to b

*xx*

### ("viscosity"

### dissipation), then the wave equation describes a train of solitons with different am-

### plitudes and this gives the s t r u c t u r e of a collisionless shock [ *Sagdeev 1966]. T h e next *

### section is devoted to the analysis of the wave equation, taking the effects of weak

### mass loading into account, t h a t is, the above mentioned ratio of spatial scales a p / i ? £ ,

### is small. We will see that the main terms have exactly the form of "hydrodynamic"

### as well as of "viscosity" dissipation.

**5. Q U A S I S T A T I O N A R Y P L A S M A M O T I O N S : ** **E V O L U T I O N OF N O N L I N E A R WAVES **

### Now we consider the slow evolution of quasistationary nonlinear waves in a hot plasma flow with mass loading, that is we are interested in the nonlinear solution of the set (2.3-4),(2.8-11) which is quasistationary in the reference frame of the neutral gas. In the absence of mass loading such a solution has the form of a soliton, as we saw in the preceeding section, with velocity u in a frame of reference moving with the plasma and with velocity f — u relative to the cometary nucleus. If one takes weak mass loading into account we must add to (4.1) dissipative terms which can be determined from the set (2.3-4),(2.8-11). The main term, in the case where the ratio of dispersive scale to mass loading scale is small, is proportional to b and to b

^{xx}*. We thus find the nonlinear equation *

### o / t/2 \

*b**x** + 3bb**x** + a**2D**b**xxx** + u**v**b**2 *

*ii \*

### \ / u \ I 6 = 0,

*p(f-u)*[\\fj f J rnuy\fj V ƒ, *

### (5.1) where the viscosity dissipation term u

^{v}

### has a rather complicated dependence on the large scale motion.

### In the case of weak nonlinearity

(*b / a o b*

*x x*

### <C 1) and weak mass loading, (5.1) can be solved by perturbation techniques (adiabatic motion of the soliton (4.3) with slowly varying amplitude). As customary for a perturbation method, the con- dition that the zeroth-order solution (4.3) be orthogonal to the first-order solution determines the evolution of the soliton amplitude:

2

### 1.(0) , * um* *f *

*—*

### o} H

*f-u *

### From this equation it is evident that nonlinear waves can propagate quasistationary in a mass loaded plasma flow, provided the term between curly brackets vanishes in (5.2). This condition determines a second relation, in addition to (4.4), between the soliton amplitude and its propagation velocity. Both together give the quasistationary soliton amplitude at any point in the mass loaded plasma flow. Using the large scale motion characteristics (3.2), we can express the only parameter in the dissipative term of (5.2), which depends on these parameters, in terms of the magnetosonic **Mach number of the plasma flow M = ƒ / V : **

**Mach number of the plasma flow M = ƒ / V :**

*ht* ^{=} * W*

^{=}

^{ ( 5 3}

^{) }

### mu 2(1 - M

^{2}

### ) '

^{ {}

### '

### One then finds that the dissipation term is proportional to 1 / ' / 5

### 1 - M

^{2 }

**where £ = 1 — u/f. This second-order polynomial in f has two roots **

**where £ = 1 — u/f. This second-order polynomial in f has two roots**

### 6 -

**M**

**M**

^{l}** - 1 ** **Z2**

**- 1**

**Z2**

^{ -}

### 2 , (5-5) **M**

**M**

^{2}** + 0.5' **

**+ 0.5'**

### that determine the propagation velocities of quasistationary localized solutions. Here we consider not only soliton-type solutions of the wave equation (5.2) but also a wider class of quasistationary solutions as a wave packet or a train of solitons, which is the limiting form of a moving collisionless shock in the case of small viscosity.

### All such solutions are the more stationary in a mass loaded plasma flow, the closer

**their propagation velocities are to the values (5.5). For linear magnetosonic waves **

**their propagation velocities are to the values (5.5). For linear magnetosonic waves**

**(£ = ± l / A f ) one can reproduce from (5.5) the result of Galeev and, Khabibrakhmanov **

**(£ = ± l / A f ) one can reproduce from (5.5) the result of Galeev and, Khabibrakhmanov**

### [1989] that magnetosonic waves propagating forward through a mass loaded plasma

**become unstable as the local magnetosonic Mach number decreases to M — 2. For **

**become unstable as the local magnetosonic Mach number decreases to M — 2. For**

**backward propagating linear magnetosonic waves this value is a little less, M — 1.74. **

**backward propagating linear magnetosonic waves this value is a little less, M — 1.74.**

**C O N C L U S I O N S **

### We have shown in this paper that, although supermagnetosonic mass loaded plasma flows with M > 2 are stable against infinitesimal disturbances, sufficiently- strong nonlinear disturbances have a definite growth rate due to mass loading. There are solutions for the wave equation of the form of l/cosh

^{2}

### -solitons with a threshold amplitude, which are quasistationary — their growth rates are zero — on the ion- ization time scale. It seems that such solitons are a permanent feature of the mass loaded plasma flow, because arbitrary disturbances of the KdV-equation split into a train of solitons, as is well known from soliton theory. According to our results, soli- tons with amplitudes below the threshold will be damped, whereas stronger ones will be accelerated and their amplitudes will grow until overturning. Only solitons with threshold amplitudes can be quasistationary. At the point in the plasma flow where the magnetoaonic Mach number equals 2, the threshold amplitude for forward propa- gating waves decreases to zero. As was pointed out by Galeev and Khabibrakhmanov [1989], this can be indicative of the beginning of the bowshock formation process in the neighborhood of comets. Similar conclusions can be drawn for the other astro- physical occurrences of mass loaded plasma flow, of which some indicative examples were given in the Introduction. These are, however, as yet less amenable to in situ observational verification.

**A C K N O W L E D G M E N T S **

### One of authors (I.Kh.Kh) is grateful to A. A. Galeev for fruitful discussions

### of many aspects concerning the present paper. F.V. thanks the National Fund for

### Scientific Research (Belgium) for a research grant. He is also grateful to M. Ackerman,

### J. Lemaire and M. Roth (Belgian Institute for Space Aeronomy) for their help which

### made the cooperation between both authors so enjoyable.

**A P P E N D I X A: D E S C R I P T I O N OF N E W L Y I O N I Z E D P A R T I C L E S ** Let us suppose that the distribution function of the neutral particles in a cometocentric frame is maxwellian with a characteristic temperature m A / 2 , m being the mass of the particles:

### We change variables to the magnetic moment fj. = Uj_/2B, with v\ — v

^{2}* + v*

^{2}*, and * to the gyrophase <j> of the ionized particles:

*v**x* = \J2\JLB c o s <J) + ƒ,

### _ (A. 2) *v*

*y*

### = \/2(j.B s'mcf) — g,

### where ƒ = E

^{y}*jB and —y — -E*

^{x}*[ B are the components of the drift velocity of the * plasma. One then gets the initial distribution function of the newly ionized particles, in the following form:

*1 ( 2fj.B + g*

^{2}* + f*

^{2}*\ ( 2V2^B ~ . - \ *

### $

0### = — exp I — 1 exp I ——(ƒ cos0 - gs\n<p) 1 . (A.3) We now see that the description of the newly ionized particles is simpler if one replaces

### the gyrophase (j> by a modified gyrophase <j> = <j> + n + a, with a = arctan(g/f) the initial gyrophase of the particle at the moment of ionization. Taking into account

### the following series expansion exp(2cos#) = -M

2### ) + 2 X^^li Ik[

z### ) cos k8, where Ik{

z### ) are the modified Bessel functions of the first kind, and using the expansions for large arguments (z » l), namely I„(z) = ^ = ( l -

^{ 4t}

### '

^{8}

### "~

^{1}

### H j , the distribution function for the newly ionized ions can be expressed in the small temperature limit in the form:

### lim — 7 = = = exp f - ^ — ( I+ 2 jr cos k<j>)

**A — » 0**

* TTy/iTrAwo \*

^{ A /}* ^ J *

## ( ^{oo }

### 1 + 2 ^ cos k<p *k=i *

### where v = y/2^B is the initial velocity of the particles perpendicular to the magnetic

### field and v

^{0}* = \JP + g*

^{2}### is the local drift velocity which the newly ionized particles

### A P P E N D I X B : V L A S O V E Q U A T I O N F O R T H E D I S T R I B U T I O N F U N C T I O N I N V A R I A B L E S (n, <f>)

### Now we consider the most simple geometry of a moving plasma. The mag- netic field vector has only one component, along the z-axis, so that B = ( 0 , 0 , 5 ) . The electric field E = ( E

^{x}

### , E

^{y}

### , 0 ) only has components perpendicular to B . In this geometry the drift velocity of the plasma will have two components: V d = {E

^{y}

### /B,-E

^{x}

### /B,0) = (ƒ,-<?, 0).

### The Vlasov equation for the distribution function F(t,r,v), describing the motion of particles in the above defined field configuration, turns out to be more tractable if one uses instead of v the variables (J. (magnetic moment) and <f> (modified gyrophase as defined in Appendix A) and expresses F throughout as a function F[t,r,fj,,<j>), so that

### dF_(W dr dF_ d»d£ d$_c)F_

### ~dt ~ dt

^{ +}

### dt ' dr

^{ +}

### dt d\i

^{ +}

### dt d<j> '

### To calculate the derivatives d^/dt and d(f>/dt one uses the particle equations of motion d

2### x _ eB dy e ^

### dt

^{2}

### m dt m d

^{2}

### y eB dx e

### dt

^{2}

### m dt m

^{ y}

### ' With the following rule for the change of variables:

### d

x### — = - yJl\xB cos{<(> - a) + ƒ,

dt

### (B.3)

### (it!

### — = -y/2fj,Bsm(<f> - a) - g, at

### where a = arctan(g/f) is the initial gyrophase of a particle (significant only for newly added particles) and using the chain rule for the derivative d/dt\

### d • d . d . d

### 37 =

=### +

+### dt dt dx dy

^{ tX}}

### ,

### d d d ( d d\ ^ '

### = dt+

^{f}

### Tx-

^{g}

### d~y-^ H* ^{ ~}

^{a)}

^{Tx}

^{+}

^{ Sin(} * ^{ -}

^{a)}

^{d-y)> }

### one obtains expressions for ji and <f>\

*ii=^(-TB + B{g*

*y*

*- f*

*T*

*)) *

### - [

c o s### ( ^ _

a### ) - f ƒ ) + sin(0 - a) + Tg) ] (B.5) + n [cos 2(0 - a) ( ~ f

*x*

* - g*

*y*

*) + sin 2(0 - a) (g*

*x*

* - f*

*y*

*)\ , *

*4> =ot - ft + ^ { f*

*y*

* + g*

*x*

*) - (cos{<f> - a)Tg + sin(0 - a ) f ƒ ) *

*+ ^{cos2{<j>-a){g*

*x*

*-f*

*y*

*)+sm2{<f>-a){g*

*y*

* + f*

*x*

*) \ . (B.6) *

### Here ft = eB/m is the gyrofrequency of a particle, and s t a n d a r d notations for the derivatives d./dx = .

*x*

### a n d d./dy = .

*y*

### are used. T h e operator T is defined as T = *d/dt + f d / d x — gd/dy. We note t h a t the first t e r m in the expression for (i vanishes * due to Maxwell's equation ? x E + d B jdl — 0, rewritten in our notationc ae

*TB = B ( g*

*y*

*- f*

*x*

### ) . (B.7)

### This simply means t h a t /x is an adiabatic invariant of the motion of the particle, a n d this is a particular incentive to use the magnetic moment in the following expansions of the distribution function for slow motions.

### T h e Vlasov equation for the distribution function F(fj.,4>), neglecting the motion along the magnetic field lines, looks like

*LF = St{F), (B.8) *

### the operator L being defined as:

*- - d 11 , , , fTg - gTf \ d *

*a fiB*

*x*

* + Tf d * *(* *i* *s* *. f g r - g f A a *

### f 5 1 5 1 +

M**c o s** 2 ( ^ - a ) | ( - /

s### - g

*y*

*)— + ~(g*

*x*

* - ƒ „ ) - | *

### / <3 1 1 + Ai sin 2(<f> - a) | ( - f

*y*

* + g*

*x*

*) — + - + /*) ^ } • *

*ƒ?«, - g fy \ d * *ƒ*

*2*

* + g*

*2*

* J d<t> *

### (B.9)

**A P P E N D I X C: E X P A N S I O N OF D I S T R I B U T I O N F U N C T I O N FOR ** **SLOW M O T I O N S **

### To solve equation (B.8) we suppose that terms in T are small compared to those in fl. T h e righthand side of equation (B.8), according to (A.4) and taking the Jacobian of the transformation (B.3) into account, will have the form:

*St(F) = - tj) + , (C.l) *

### where u is the production rate of the newly ionized particles and 77 = (ƒ

2### + g

*2*

*)/2B * is the initial value of their magnetic moment.

### T h e zeroth-order approximation of equation (B.8) gives as only result t h a t the zeroth-order distribution function Fo is independent of Averaging the equation for Fo over the gyrophase 4> according to

### 1 C

*2*

** *

### ± < . . . > « # , (C.2) 2?r Jo

### we can determine the equation governing Fo:

### T F

0### = ^ 6 { n - v ) . (C.3)

**For the next approximation we split the distribution function into two parts according ** **to F\ = Fi of which the first part is independent of <f>. Averaging now the first-** **order equation **

**to F\ = Fi of which the first part is independent of <f>. Averaging now the first-**

**f**

^{ A}* - n* **^ - V ^ B c o s - . ) { f**

**^ - V ^ B c o s - . ) { f**

^{ +}** I - f f ) **

**I - f f )**

**aFo dF**

^{0}** (**

^{C - 4}**) **

**- Aicos2(0 - ot)(/**

^{x}** + + fj-s'm2{(f) - a ) [ - f**

**+ + fj-s'm2{(f) - a ) [ - f**

^{y}* + fx)-^-*

**=27\F**

^{0}** X ]** **oo **

^{ c o s }**fc=i **

**over one obtains an equation for the part F\ which is independent of <f>: **

**over one obtains an equation for the part F\ which is independent of <f>:**

*TFi = 0. (C.5) * **Substituting the solution of this equation back into (C.4), we get an equation for F\, **

**Substituting the solution of this equation back into (C.4), we get an equation for F\,**

**with solution **

*(dF*

^{0}* 1 / „ ~ dF*

^{0 }*nF*

^{1}* = - y/2sin{4> - a) j + i [~nB*

^{x}* - f f ) ^ *

**+ V ^ B cos(* - a) { ^ + I { - , B**

**+ V ^ B cos(* - a) { ^ + I { - , B**

^{y}* + f „ ) ^ *

**(C.6) ** **- - s i n 2 ( 0 - a ) ( /**

^{z}** + 0 » ) " ^ **

**a , w , ^ ~ „ sin kd> **

**- | cos2(<£ - «)(-/„ + - 2TF**

**- | cos2(<£ - «)(-/„ + - 2TF**

^{0}* £ *

*^ k= 1 *

**The second approximation to the distribution function can found from the equation **

*n d*

*J l*

^{ =}* f h + f*

^{ F l}* - V ^ B cos(0 - a)*

^{ a ( i}**y ^**

^{2 ) }**^ . (C.7) ** *- V W s m ( < f > - « )*

*d { F l d + y F 2 )*

*. *

**Again averaging this over the gyrophase, we obtain the «^-independent part of the ** **second approximation to the distribution function: **

**nF**

^{2}** = n{g**

^{x}** + - ^/2^B ( s i n a ^ + c o s a ^ ) . (C.8) **

### Substitution of this result back into (C.7) allows one to determine the remaining part:

### n

2### F

2* = *

*dTF*

^{0}* ( nTB*

^{x}*+T*

^{2}*f n. .*

^{ nr}* A dF*

^{0 }*cos(0 - a) \J2iiB —^ +* *dx \ B*

^{ x}*—* 4

^{ J}*-.+ Z (4g*

^{ v}

### " "

^{xy}

^{ Jyy}* - ƒ « + 3 f*

^{y y}*' J dii * *) *

### ^2 .

*^ /—~ dTFo ƒ iiTBy — T g p . .*

^{ n}* , , \ dF*

^{0}*\ * *+ sin(0-a)V^Bj-^+ I -Z*

^{ y}*-*

^{ y}*- + ^(g*

^{yy}*-Zg*

^{xx}*-4f*

^{xy}*) J \ *

+

### f c o s

2 (### , - -

+### 1 ( „ < * . - B „ ,

+### § f < *

+### , . , ) f }

### i * c) W *

*+ \/2fiB — ~ ~~ cos 3(0 - a) - 2<7.*

^{T!/}

### - ,/

^{Tr}

### ) 12 dfj.

### 3 ri

### - Sin 3(0 - a) (g

^{yy}

### + 2 /

^{I y}

### - - 2 M S c o s ^ - a ) | c o s a —

+ S i n ( i### — |

*„ , , * f - d*

^{2}*Fo d*

^{2}*Fo \ * + 2//Bsin(0 - a) j s i n a - ^ - +

c o s a### d

x### d y J

### „ cos * kd> *

*+ 2T*

^{2}

*FoY^ *

*k=l *

k 2
*.k =*

^{ 2 }### oo

### . — r j r ^ z cos((A; - 1)0 + a) cos((fc + 1)0 - a ) \ cos(20 - a) [ d T F

0### - ** \** ^{ ^ ( +}

k { k + 1 )** ) + > **

### , — - J ^ /

S*in((fc — 1)0 + a ) sin((fe + 1)0 — <x)\ sin(20 - a) [ dTF*

*0*+

### ^

B### { I j I k(k - 1 ) w + 1 ) ) 2 ƒ ~ W

### (C.9)

### From the equation for the third-order approximation to the distribution function:

**nd**

**J±=TF**

**J±=TF**

**z**** + TF**

**+ TF**

**2**** - ^Bcos{4> - a) _ yj2\iB sin(<£ - a)**

**- ^Bcos{4> - a) _ yj2\iB sin(<£ - a)**

**d****SEl±IA ^ **

**SEl±IA ^**

**o<p ox ay ** (C.10)

**o<p ox ay**

### we need only the part averaged over <f>:

+

### j ^ c o s a - ^ s i n a j ( C . l l )

**HB f fd**

**HB f fd**

**2****F**

**F**

**0**** d**

**d**

**2****F**

**F**

**0****\ d*F**

**\ d*F**

**0**** . . \ ** cos 2a — 2 sin 2a > ,

**. . \**

### 2 I \ dx

**2**** dy**

**dy**

**2**** J dxdy **

**J dxdy**

### and the first harmonics in the gyrophase <f> of the (/»-dependent part:

O , ,—=\dT2Fo 3ßB ( d3F0 d3F0\

n F3 = r + mt)

( f3F0 Sy/2JIB f d2TF0 d2TF0

-f cos a 2—== . 1 r—T

\ yfïixB 4 y dx2 dy2

3 / . df2F0 o df2F0V + - sm2a— cos 2a—

2\ dy dx J

UB ( ' (d3F0 d3F0 \ . d3F0 \

+ T r s 2 a { - d ^ - a * * ) -2sm2ad^d~y)
1 / f f d2TF0 d2TF0\ n . „ d2f \
+ - V W [ c o s Z a ^ - ^ - - ^ - J - 2^{S 1}n 3 a — J

/ f3/ b ^ , ÔFû

+ ( ^ s s ^ B * - ^ - + Y{ B x y y + B x x x ) + T T { f y y + 9 x y ))

U \ '/m\df2Fo Z ß B ( d Z p° 4. 5 3 M - cos(<£ - a ) v / 2 m 5 I — + ^ r j

/ f3F0 3y/2jïB ( d2TF0 d2TF0

-sm a [27^Ë - [ +

3 f n dT2F0 . 0 df2F0\

+ ï [ c o s 2 a - ^ r + s m 2 a - d ^ ) .

uB / ( d3F0 d3F0\ . d3F0 + T [COS2a { d ^ d - y - ^ ) - 2 S m 2 a d ^

1 > f ( d2TF0 d2TF0\ n 0 d2f \ + - ^ s i n 3 a + J - 2 c o s 3 a — j

+ + f-W + T + * » " > - T f + m ) f ( C . 1 2 )