Plasma expansion into vacuum — A hydrodynamic approach

Plasma expansion into vacuum — A hydrodynamic approach

PLASMA EXPANSION INTO VACUUM A HYDRODYNAMIC APPROACH - Ch. SACK and H. SCHAMEL Institut für Theoretische Physik, Ruhr-Universikit Bochum, D-4630 Boch...

6MB Sizes 0 Downloads 47 Views

PLASMA EXPANSION INTO VACUUM A HYDRODYNAMIC APPROACH -

Ch. SACK and H. SCHAMEL Institut für Theoretische Physik, Ruhr-Universikit Bochum, D-4630 Bochum 1, West Germany

I

NORTH-HOLLAND AMSTERDAM -

PHYSICS REPORTS (Review Section of Physics Letters) 156, No. 6 (1987) 311—395. North-Holland, Amsterdam

PLASMA EXPANSION INTO VACUUM

-

A HYDRODYNAMIC APPROACH

Ch. SACKt and H. SCHAMELtt Institut für Theoretische Physik, Ruhr-Universitdt Bochum, D-4630 Bochum 1, West Germany Received June 1987

Contents: 1. Introduction 2. The plasma expansion model 2.1. Basic equations, initial- and boundary conditions 2.2. Constraints on the electron thermodynamics 3. Plasma expansion within the framework of gas dynamics 3.1. Self-similar theory 3.2. Simple waves and Riemann invariants 4. Numerical methods and analysis 4.1. Lagrangian equations of ion motion 4.2. Two iterative methods for the nonlinear Poisson’s equation 4.3. Boundary conditions 4.4. Conservation laws 4.5. Variational principle 5. Numerical results for the dissipationless plasma expansion 6. Theory of bunching and wave breaking in ion dynamics 6.1. Scalar wave equation 6.2. Perturbative solution of the scalar wave equation for isothermal electrons

313 316 316 319 322 322 327 332 332 334 336 339 342 343 351 352

7. Comparison with previous hydrodynamic and kinetic models 8. Navier—Stokes viscosity and implicit Lagrangian scheme 9. The early stage of viscid plasma expansion 10. The long-time behaviour 10.1. The debunching process 10.2. The asymptotic velocity of fast ions 10.3. Intermediate asymptotics 11. Comparison with experiments 12. Summary and conclusions Appendix A: Bunching and wave breaking in nonlinear Langmuir oscillations Appendix B: Spatial discretization of Poisson’s equation Appendix C: Implicit difference scheme for the equations of viscid ion motion Appendix D: Conventional stability analysis of the discretized set of equations References

360 364 367 371 371 373 374 377 380 381 385 388 390 393

354

Abstract: Based on the hydrodynamic description, the planar expansion of a plasma into vacuum is investigated numerically and analytically. Several dynamical structures are found and explained. In the inviscid case, ion wave collapse is the most striking feature. Viscosity prevents the collapse and allows long-term calculations. Three phases in the evolution can be distinguished. In the time asymptotic regime the self-similar state is approached. This approach is preceded by a state of intermediate asymptotics and applies to charge separation as well.

Now at JET Joint Undertaking, Culham Laboratory, Abingdon, Oxon 0X14 3EA, England. ~ Now at Physikalisches Institut, Universität Bayreuth, D-8580 Bayreuth, Postfach 3008, West Germany.

Single orders for thi.s issue PHYSICS REPORTS (Review Section of Physics Letters) 156, No. 6 (1987) 311—395. Copies of this issue may be obtained at the price given below. All orders should be sent directly to the Publisher. Orders must be accompanied by check. Single issue price Dfl. 61.00, postage included.

0 370-1573/87 / $29.75

© Elsevier

Science Publishers B.V. (North-Holland Physics Publishing Division)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

313

1. Introduction The investigation of a spontaneously released or produced plasma expanding into vacuum is of considerable interest to the general understanding of unsteady, nonlinear behaviour of ionized gases. The special role plasma expansion is playing as a basic process in astrophysical flow problems has been repeatedly emphasized [1, 2]. In fusion research expansion processes are just as fundamental if one thinks of phenomena such as the arcing on surfaces of limiters and divertors [31, the ablation of solid hydrogen pellets in refuelling of tokamaks [4, 5], and the laser—pellet interaction in inertial confinement. In the latter case, e.g., a detailed knowledge of the dynamics of the coronal plasma produced by the laser beam is required in order to estimate its effects on the implosion efficiency [6—12].As it is well known, the existence of a component of high energetic ions deteriorates the pellet compression. Although these fast ions represent only a relatively small fraction of the total ablated mass, they carry off up to 50 per cent of the absorbed laser energy [8, 131. The feasibility of laser-induced fusion, therefore, is essentially coupled to the minimization of the number of fast ions which presumes a basic understanding of their generation mechanism. Accelerated ions in connection with expanding plasmas are observed in other areas, too, such as in the polar wind [14, 151, in cathode flares [161,in vacuum arcs [17] as well as in exploding wires [18, 19]. As early as 1930 Tanberg [201found plasma jets of high speed in a pulsed gas discharge being ejected from the plasma region in front of the cathode. An explanation for this effect based on the mechanism of ambipolar ion acceleration by electrons was not given until the early sixties by Plyutto [211,and by Hendel and Reboul [22]. Due to the higher mobility of the electrons, an electric (ambipolar) field is generated which prevents the electrons from escaping and, at the same time, accelerates the ions towards the vacuum or the less dense medium. Experimental [23—3 1] and theoretical [32—521 investigations carried out during the following years have further contributed to a clarification of the ion acceleration process, and have discovered nonlinear phenomena being involved in the plasma expansion such as shock waves, wavebreaking and anomalous dissipative effects. According to the general plasma theory, expanding plasmas are described by both kinetic and hydrodynamic models. Completed by Maxwell’s equations, quite complicated sets of equations are obtained which cannot be solved analytically. Even the numerical solution of these sets causes a lot of difficulties. However, in order to comprehend basic processes in plasma expansion, it is often sufficient to consider simplified models being suggested by experimental observations. Such a simplified model is the one-dimensional planar expansion of a rarefied magnetic field-free plasma into vacuum. The space-time evolution of the rarefied plasma is described by a system of collisionless kinetic equations for electrons and ions in connection with Poisson’s equation for the electrostatic potential. Assuming that the electrons always stay in equilibrium with the electrostatic potential and that the plasma behaves quasi-neutral, the set of equations is reduced to one single equation, the Ion—Vlasov equation which is, however, nonlinear. Because of the quasi-neutrality assumption this equation can be transformed to self-similar variables with reference to usual hydrodynamics (see, e.g. [53,54]). This means that in the solution of the problem the time t and the space-coordinate x appear only in the combination x/t. In the case of Boltzmann-distributed isothermal electrons such solutions of the plasma expansion problem were firstly presented by Gurevich et al. [33]. They indicate a strong, continuous ion acceleration driven by the ambipolar electric field which has been found later on in other extended investigations, too [37, 38, 42, 46, 49, 50]. Coincident with the acceleration the effective ion temperature rapidly cools down, and finally the motion of the accelerated ions can be described by the hydro-

314

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

dynamic equations of an ideal, isothermal gas [33]. Assuming cold ions from the outset, the self-similar solution of the hydrodynamic equations yields at a fixed time t a density profile exponentially decreasing, and a velocity profile linearly increasing in x [36, 43]. An unphysical feature of this solution is that the electric field is singular at t = 0 and does not fall off for x + ~ at any other time. Because of these shortcomings and because of the fact that a characteristic scale length is missing, the self-similar theory is only conditionally suitable to describe the plasma expansion for an arbitrary initial- and boundary value problem. In which sense the quasi-neutral self-similar solution is approached in an expanding plasma, will be studied among others in this review-kind article. A more general description of the nonlinear motion of a collisionless quasi-neutral plasma, which includes the self-similar solutions as a special case, is based on the Riemann solutions for so-called simple waves of arbitrary amplitude in usual hydrodynamics (Riemann (1860), compare also [53]).The most noteworthy property of this solution is that any compressional wave progressively steepens, and finally breaks after a finite time. In this case the flow profile becomes multiple-valued in space, and the equations loose their global validity [53]. Whereas in usual hydrodynamics the occurrence of singularities in the flow profile is prevented by dissipation, dispersive effects come into play in the plasma case. They are believed to limit the steepening of finite amplitude waves and inhomogeneities, leading to the development of solitary and/or periodic waves [55]. Dispersive effects result from the deviation from quasi-neutrality [32, 34, 35]. In this case the Debye length is the characteristic scale length of the problem, and the potential has to be calculated from Poisson’s equation. By this means, the singularity of the electric field at t = 0 existing in the quasi-neutral self-similar solution is removed. There has been several attempts to simplify the system of hydrodynamic ion equations and Poisson’s equation by appropriate self-similar transformations and thus, to render them accessible to an at least approximate analytical treatment [44, 48]. However, the general initial- and boundary value problem for an expanding plasma including charge separation can only be solved numerically. For the numerical investigation of plasma flow on a hydrodynamic basis finite difference schemes, such as, e.g., Eulerian- and Lagrangian schemes, have been developed. In this context repeated interest has been devoted to the isothermal expansion of a cold ion plasma into vacuum [37, 42, 49, 50]. Besides a lot of according results on this subject being confirmed by particle simulations of the problem [46], and being discussed in detail in this review, quite contradictory statements have been recorded with respect to the concrete space-time behaviour of the fastest ions. Especially it was not clear whether the occurrence of a density hump at the ion front shown, e.g., by Crow et al. [42], and by Gurevich et al. [49], is due to a numerical effect or whether there is a physical phenomenon hidden behind it. True et al. [50],who solved the same set of equations, simply state they had not observed such a density hump. In an earlier paper of the authors [56], the one-dimensional planar expansion of a collisionless plasma was studied under the influence of an obliquely incident electromagnetic wave. Using an explicit Lagrangian scheme for the ion hydrodynamics, a local, explosive increase was found in the ion density leading to a sudden breakdown of the numerical solution. A similar breakdown was also seen when the radiation field was switched off. In subsequent papers [57, 58], the present authors investigated this pure plasma expansion problem in more detail. It was found that the density hump and the breakdown of the hydrodynamic solution are of the same origin being connected with the appearance of a simple wave structure in the expansion problem which develops from an extremely inhomogeneous initial density. The main topic of this review is to evaluate the nonlinear structures arising in the plasma expansion modelled by the hydrodynamic equations of motion. The paper is organized as follows: —~

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

315

In section 2 the mathematical formulation of the problem and its physical foundation is presented. Section 2.1 contains the basic equations including the general initial- and boundary conditions. In section 2.2 it is shown that in the case of a polytropic electron equation of state the specification of the boundary conditions in Poisson’s equation imposes constraints on the electron thermodynamic behaviour. The analogy of the plasma expansion model to the usual hydrodynamics of an ideal neutral gas is discussed in section 3. Section 3.1 deals with the analytical treatment of the plasma expansion within the framework of the self-similar theory. Section 3.2 contains, the more general formulation of the expansion problem by means of the theory of simple waves of finite amplitude and of Riemann invariants. In section 4 we present the foundations of the numerical solution methods for the dissipationless hydrodynamic equations and for Poisson’s equation as well as some corresponding analytical peculiarities. Section 4.1 contains the transformation of the Eulerian equations of motion to the Lagrangian picture, whereas section 4.2 is devoted to the iterative solution methods for the nonlinear Poisson’s equation. In section 4.3 the boundary conditions for the ion equations and for Poisson’s equation are adapted to the numerical solution methods. The conservation laws resulting from the basic equations are deduced in section 4.4. In section 4.5 it is shown that the ion hydrodynamic equations can be completely derived from a variational principle under the constraint of Poisson’s equation. In section 5 the numerical results of the dissipationless plasma expansion are presented. The effects of charge separation and of quasi-neutrality on the global behaviour of the solution are discussed, as well as the influence of electron thermodynamics. In those cases where the solution breaks down, a sharp spike in the ion density emerges, increasing explosively. At the same time the electric field and the ion velocity steepen without limits. These properties of the solution being analytically analysed in detail in section 6, suggest the occurrence of wavebreaking in connection with density bunching as it is known for Langmuir waves and for ion waves of finite amplitude. In section 6.1, a single scalar wave equation is derived describing the whole ion dynamics including charge separation. Its solution by perturbative means in the vicinity of the breaking point given in section 6.2 reproduces all the properties found numerically. Thus, it is proven that the singularity in the ion dynamics is not a numerical effect but an expression of the extreme nonlinearity of the plasma expansion into vacuum. In section 7 comparisons are made with previous hydrodynamic and kinetic investigations. The kinetic calculations for an equivalent problem give a hint which types of phasespace structures can be expected after wavebreaking. Within the hydrodynamic description wavebreaking is prevented by including dissipative effects. For this purpose a viscosity term of the Navier—Stokes type is introduced in the ion momentum equation rendering possible the implicit formulation of the difference scheme. In section 8 the properties of such a scheme, especially its numerical stability, are discussed. The full-implicit scheme turns out to be in favour of the stability behaviour of the code allowing long-term calculations of the expanding plasma without showing any numerical instability. The results of these calculations are presented in sections 9 and 10. Section 9 deals with the early stage of viscid plasma expansion where the initial acceleration and the ion bunching occur. The ion bunching process is now stabilized by the viscosity term and manifests itself in a pronounced density hump propagating at a constant speed of several times of the ion sound velocity. The quantitative properties of the acceleration phase and of the bunching phase depend on the viscosity, on the electron thermodynamics, and on the initial conditions.

316

Ch. Sack and H. Schamel. Plasma expansion into vacuum — A hydrodynamic approach

In section 10 the long-term behaviour of the expanding plasma is discussed. Section 10.1 shows that the two phases explained in section 9 enter a third phase which is initiated by the debunching process. The debunching leads to a rapid decrease of the density hump and gives rise to a further acceleration of the ion front. As presented in section 10.2 the fast ions at the front achieve a velocity which is given by the self-similar theory. The self-similar behaviour of the plasma expansion as a state of intermediate asymptotics is exemplarily proven in section 10.3 by a perturbative analysis of the single scalar wave equation for the isothermal quasi-neutral case. In section 11 the results of our numerical calculations are used to interpret the measurements of some selected plasma expansion experiments being especially contradictory with respect to the velocity of the fastest ions. Section 12 contains the summary of this review and the conclusions. In the Appendix the nonlinear behaviour of Langmuir waves and the details of the finite difference scheme and of the stability analysis are treated separately.

2. The plasma expansion model 2.1. Basic equations, initial- and boundary conditions As mentioned in the introduction, we investigate the one-dimensional, planar expansion of a strongly inhomogeneous, magnetic field-free plasma into vacuum. It is supposed that the plasma consists of single charged ions with mass m~and charge +e, and of electrons with mass me and charge —e; e represents the elementary charge. Assuming that the ion temperature T1 is negligible compared to the electron temperature T~,e.g. T1/T~—~0, we have the following two-fluid equations: +

~i~(n1v)= 0

(2.1)

+

v~a~v~ = —-f--

(2.2)

ax(neve) = 0

(2.3)

+

e c31V~+ Ve i3xVe

=



I —



(2.4) (2.5)

The equations (2.1)—(2.4) describe the space-time development of the density n and of the fluid velocity o of the two-component plasma, where the indices i and e distinguish ion and electron quantities. Equation (2.5) for the electron pressure is obtained by inserting the equations (2.3) and (2.4) into the electron energy equation; We represents the electron heat conduction. The set of equations (2.1)—(2.5) is supplemented by Poisson’s equation for the electrostatic potential 4, which couples the ions to the electrons: =

4~re(n~ni). —

(2.6)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

The equations (2.1)—(2.5) are not yet complete because an equation for achieved by assuming a polytropic equation of state for the electrons: Pe1~T~~’~k=const.,

We

317

is missing. Closure is (2.7)

where y is the polytropic exponent. The polytropic equation of state corresponds to a nonvanishing divergence of the heat conduction We~ given by axwe=—(3—y)kaxve=—(3—y)peaxve.

(2.8)

Inserting (2.8) into (2.5), one easily gets (P 5\_ d

(Pe~

(2.9)

~

from which follows by integrating (2.9) Pe~12~ n~0

(2.10)

where Peo and n~0are the corresponding quantities of the electron pressure and of the density in the unperturbed homogeneous plasma. Assuming that the electron fluid behaves like an ideal gas, i.e. PehleKTe

(2.11)

where K is the Boltzmann constant, we get from (2.10) and (2.11) a relation between the electron temperature and the electron density: Ten~~~ const.

(2.12)

In the isothermal case, y = 1, it follows from (2.12) Te = const. as it should be. Owing to the high thermal conductivity of the plasma, the electron temperature will stay constant during the expansion process. This standard case is described by an isothermal equation of state, ‘y = 1. In certain circumstances, however, the heat transport can be reduced especially when regions of strong gradients or extremely dilute plasmas are involved. We, therefore, allow in our general description also for nonisothermal equations of states, y 1. In fact, equations of states deviating from isothermality have been identified in particle simulations (see e.g. Denavit [46]) and in laboratory experiments (see e.g. Hendel and Reboul [22] or Eselevich and Fainshtein [30]). Being interested in processes taking place on the ionic time scale, we normalize the system (2.1)—(2.4) and (2.6) as follows: a~n+ 9~(nv)= 0 +

u 9~v= —a~4

(2.13) (2.14)

318

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

atne

8x(neve) = 0

+

(2.15)

m

1

—s (a,v~+

Ve axve)

=





a~n~

(2.16) (2.17)

where t, x, o, n and

4’ are normalized by

w ~‘, AD, c~,n0 and K Teo/ e. The latter quantities are the ion 2/m 1~, the electron Debye length AD = (K Teo/4lTnoe2)”2 and the ion plasma frequency ~ Ten/mi)112, = (4~rn0e 1) the index 0 refers to the unperturbed, homogeneous plasma sound speed c~= (K where (nen = As long as the ion fluid velocity is much less than the electron thermal velocity °the = (K Teo/ me)’ / 2 noting that Vthe/CS = (m 112 1 in these units, we can neglect the left-hand side in the electron 1/m~)which then reads momentum equation (2.16), ~‘

(2.18) Integrating eq. (2.18) yields the electron density ne as a function of the electrostatic potential =

(~+7

—1

4:

l/(y-l)

(2.19)

.

From the requirement that the electron density falls off from unity in the unperturbed plasma to zero in the vacuum, it follows that the potential is restricted to the interval 0 and —‘y/(y 1). In the isothermal case, ‘y = 1, we obtain from (2.19) the well-known Boltzmann law: —

ne(4)

=

e’~

(2.20)

.

For y = 2 we get a linear relation between the electron density and the potential: ‘~e(4~)=

1

+

cb/2.

(2.21)

The set of equations which will be the basis for the following investigations, therefore, reduces to: (2.22)

a~n+a~(no)=0 +

v

=

0~4= ne(4)

—a,~ —

n

(2.24) —

=

(i

+

(2.23)

1

~ e4,

~

7 1 y=1

(2.25)

We note parenthetically that this system can be derived from a variational principle which will be given in section 4.5.

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodvnamic approach

319

The solution of the equations (2.22)—(2.25) is subject to the initial and boundary conditions: a) initial conditions (t = 0): n(x,0)=n0(x),

~(x,0)=4.0(x), v(x,0)=u0(x);

b) boundary conditions (x

—~

(2.26)

±so):

v(—co,t)=0

(2.27)

7

n(+cc,t)=0,

~(+so,t)=— 1

,

v(+co,t)=0.

The potential 4~i0(x)consistent with the initial ion density n0(x) is obtained by solving Poisson’s equation (2.24) at t = 0. In contrast to Widner et al. [37], Crow et al. [42] and Gurevich et al. [49], who have used a discrete density step for the ions, we prefer a density profile with a finite transition length 1: n0(x)

=

~ arctg[exp— (x7v0)],

n0(x0)=0.5;

(2.28)

for x—~±so, n0(x) reaches the values given in (2.27). A diffuse density profile corresponds to a more realistic initial situation than a density step, having in mind, e.g., the laser—pellet interaction. In this case, the laser prepulse produces a finite density gradient through heating, ablation and ionization. The parameter 1 in (2.28) allows to control the steepness and the transition width of the density profile; for l—~0 (2.28) degenerates to a density step. In our calculations we have chosen 1 = 4 and 1 = 1, respectively. This choice, especially 1 = 1, guarantees a sharp density transition, in front of which the ion density is negligible. In most of the cases treated in this review, the ions are assumed to rest initially, i.e. v0(x) 0. In section 9 also two cases with v0(x) ~ 0 are considered. Mathematically the equations (2.22)—(2.25) are a coupled nonlinear system of partial differential equations, the general solution of which can only be obtained by numerical means. In special cases, however, one can find analytical solutions of the system or one can at least study some properties of possible solutions analytically. In section 3 we shall come back to this point. In the following subsection we present a special analytical solution of Poisson’s equation and discuss the consequences on the thermodynamic behaviour of the electrons. 2.2. Constraints on the electron thermodynamics From earlier numerical investigations of the plasma expansion into vacuum taking into account charge separation (see, e.g. [37, 42, 46, 49, 50]), it is known that the ions acquire a well-defined front, beyond which there exists a pure electron cloud. We may suppose the same in our case, and introduce x~,the time-dependent position of the ion front. In this case Poisson’s equation is given in the electron cloud region by =

=

(i



+

~

1

1/(y-l)

where we assume for the moment ‘y

x > x~,

,

1.

(2.29)

320

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

Multiplying eq. (2.29) with çb’ and integrating, we obtain: —1 E2 ~ =(i+~ as~,

~) ~

(2.30)

where we use the boundary conditions 4’(x—* +so) = 0 and ç5(x—~+co) = —y/(y 1). Taking the square root of (2.30) and replacing the electrostatic field E by we are able to perform a further integration: —

—~‘,

1

[(~



[x x~(t)]= —

4,

Solving (2.31) for

~

+ ~

~

1

(2~’2~~1) ~



]~(()

we get the potential within the electron cloud x > x~: 2

~(x,

t)

=



(2.31)

~

fi



A2~’~[i

—2ey—1)/(2—y)



+

~

A~2~(x xf(t))] —

with A=[1+7



1

}

(2.32)

~t(t)]

where 1 takes the value 4~at the ion front Xf. The time dependence of the potential in the electron cloud region enters through xf(t) and 4 1(t) which are analytically unknown, and which will be provided by the numerical solution treating the full problem. If one assumes a step profile for the ions Ji, n(x,0)—i~0

~~0nsi~(t=O)

one gets ~(x0, 0) 4~(t= 0) = —1. This result follows from the continuity of çb and E at the density jump x = x0 and has already been formulated by Crow et al. [42] for the isothermal case y = 1; it is, however, valid for any ‘y. From eq. (2.32) we obtain the electric field E = —4.’ and the electron density ‘~‘e= [1 + ((‘y— 1)/y) ~]1/(Y-1): E = V’~A [i

2-y)

+

A~2~(x xi)] -y/( 2ç~—~ —



2 =

A21~[i

+

ç~-~ A~2~(x

with y ~ 1 and y ~ 2. The case y ~=~t21n[i+~



=

(2.33)

—2/(2—y)

xi)]

(2.34)

1 was obtained by Crow et al. [42]:

(x_x~)]

(2.35)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

E = ~exp(~) [i+ exp(~/2)(x

~,1 Li+

nee

For ‘y

=

exp(q5~/2) ,~



12 (x—xf)]

xi)]

=

E

—2{1

(1 + çbfl2) exp[—(x



\/‘~ (1 +

=

~~/2) exp[—(x

(1 + ~~/2) exp[—(x







(2.36)

(2.37)

.

2 we get an exponential dependence of

321

~,

E and ne:

x~)!V~]}

(2.38)

x~)/\/~]

(2.39)

xf)/\/~].

(2.40)

In the following we investigate the compatibility of the equations (2.32)—(2.34) with the boundary conditions for x—~~so (see eq. (2.27)). All equations (2.32)—(2.34) contain the term r

Ii i’

~,

~



r—

V

L

~‘

7

I

2Y)17( ~

A(

‘~

X~I J

which must go to zero for x—~+so, in order to satisfy the boundary conditions 4(x—*+co)=—y/(y—i),

for all times

t

E(x—*+cc)=0

and

n~(x_~~+so)=0

(2.42)

0. This is the case when the polytropic exponent y lies in the interval

l
(2.43)

because the exponent 11(2 y) in (2.41) remains negative, and the expression (2.41) vanishes for x—* ~ In the special cases ‘y = 1 and y = 2, the solutions (2.35)—(2.37) and (2.38)—(2.40), respectively, also satisfy the boundary conditions (2.42), so that the inequality (2.43) can be extended to —



l~y~2.

(2.44)

For y >2, however, the exponent in (2.41) becomes positive, and the boundary conditions cannot be satisfied anymore, because in the equations (2.32)—(2.34) 4, E and tie diverge or become imaginary for x—> +so. This unphysical feature of the solution can be avoided by assuming an electron front at finite x. To prove this statement, we consider (2.41) for y >2: r Ii

L

~



\,/~

2 45



~,X

X~

.

The concept of the plasma expansion into vacuum, which implies the monotonic vanishing of ne and E, can only be maintained as long as the expression in the square brackets is non-negative. From 1— ~

322

Ch. Sack and H. Schamel. Plasma expansion into vacuum — A hydrodynamic approach

we obtain the position of the electron front x~1: Xet

=

Xf +

2 ~

y >2.

(2.46)

For y >2, the boundary conditions (2.42) must be modified as follows: 7—1’

E(x=xCf)=O

and

fle(X~et)O

(2.47)

We thus arrive at the conclusion that there exists a unique relation between the thermodynamic behaviour of electrons, being described by the polytropic exponent y, and the initial and boundary conditions. The restriction of y to the interval 1 s y ~ 2 implies that the electron front is at infinity already at t = 0. For y > 2, the electron front is in the finite region, which is also true if one uses an electron equation of state describing more sophisticated trapped electrons (see [59]). The numerical methods used in this review require the restriction of the polytropic exponent to the interval described by inequality (2.44) (see section 4.3). Some results in this subsection, especially eq. (2.30), turn out to be useful in establishing the boundary condition for the numerical scheme as also shown in section 4.3. In the next section we investigate the system (2.22)—(2.25) with respect to its analytical solubility in the quasi-neutral regime. The framework of ordinary hydrodynamics will provide a convenient basis for the derivation of some properties being recovered in the numerical solution.

3. Plasma expansion within the framework of gas dynamics 3.1. Self-similar theory

The formulation of the plasma expansion on a hydrodynamic basis suggests a relation to ordinary gas dynamics. This relation is most evident if one first assumes quasi-neutrality, —

tie=~=(i+~

1

~)

l/(y-i)

and, thus, expresses the potential 3,n

+

a~(nv)= 0

(3.1)

,

4

in (2.23) in terms of the ion density n: (3.2)

2a~n. (3.3) o = —yn~ The equations (3.2) and (3.3) are identical with those describing the dynamics of an ideal neutral gas. In order to solve these equations one may, therefore, apply the known methods derived for Eulerian hydrodynamics, such as dimensional analysis or similarity treatment, which can be found, e.g., in Landau—Lifschitz’s text book [53]. However, we are not solving the equations (3.2) and (3.3) by a dimensional analysis, but shall use a general mathematical formalism based on the theory of transformation groups for partial differential equations (Ames [60], Lonngren [61]). Within this theory it is possible to derive self-similar variables and to investigate the self-similar behaviour of partial differen+

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

323

tial equations. The independent variables x and t are combined to a new variable reducing the set of partial differential equations to an ordinary one. The solutions of the latter are called self-similar, because by rescaling the solution in x-space the solution for different times can be mapped onto itself. The connection of the similarity transformation to the theory of Lie groups is given as follows (see Lonngren [61]): Similarity variables are identical with the invariants of a certain one parameter or more parameter group of transformations. The most general group is called “infinitesimal group”; it contains all possible similarity variables. The derivation of similarity variables from the infinitesimal group, however, is algebraically very cumbersome (Lonngren [61]). We restrict ourselves henceforth to the definition of a one parameter linear group G:

I t=aa~t, G:=~_ 1n, ~n=a

v=a2v

~

-

(3.4)

,

where a is a positive, real constant. The exponents a, and f3~~ i = 1, 2, have to be determined such that the equations (3.2) and (3.3) are “constant conformally invariant” under the group G [60,61]. Definition: A function F(y) is said to be “constant conformally invariant” (CCI) under G, when F(y) =f(a) F(~),

(3.5)

where f(a) is any function of the parameter a. If f(a) conformally invariant” (ACCI).

1, F(y) is said to be “absolute constant

We insert (3.4) and the transformed derivatives =

aa,

a

=

a~3~

(3.6)

into the hydrodynamic equations (3.2) and (3.3): aa1~,

a~+ aa2~,~2a~(n~Y) =0

a~1I323~s~ + a2~’~2~2

(3.7) 2)~2Pi

1j ~j

y ñ~’~ a~.

(3.8)

_ai(~

In order that the equations (3.7) and (3.8) are CCI (inclusively ACCI) according to the definition given above, it must hold: (3.9) a 1

f~2

a2 —2/32

a2—(y— 1)13k.

(3.10)

The nontrivial solutions of the equations (3.9) and (3.10) generate the similarity variables which are at the same time invariants of the group G [60,61]. Firstly we determine these invanants and eliminate the time variable t. From the first equation in

324

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

(3.4) it follows: a

and

= (~/t)’~’

x

(1/t)a2/~1x,

a

=

const

(3.11)

.

with a1 0. For a, = 0, the case we are not considering here furthermore, we obtain another set of invariants (Ames [60]). With (3.11) we obtain the first invariant =

x/ta2~1 r

=

r(x,

(3.12)

t).

In a similar way the second set of equations in (3.4) yields the further invariants n(x,

t) /tP1~a1 =:

t~(r),

(3.13)

o(x, t) /~P2~1 =: i~(r).

From the equations (3.9) and (3.10) we get the following relations for the exponents a2/a1,

/31/a1

and

f321a1: ~21a1a2~11~~1, f31/a,

=

2(A



(3.14)

A—a2/a,

1)/(y —1),

(3.15)

where y 1 is assumed. In order to determine these exponents uniquely, we require that the boundary conditions of the expansion problem formulated originally in x and t, are transformed without change onto the T-dependent variables. For x —so the boundary condition for the ion density reads —~

(3.16)

n(x—*—so,t)=1.

Because n(x—~—so, t) transforms through (3.13) to i~(r—*—so), (3.16) demands, that /31/a1 = 0. From (3.14) and (3.15) we thus get f32/a1 = 0 and A = 1. With these results the self-similar variables are obtained from (3.12) and (3.13) as follows rx/t,

n(x,t)=n(r),

(3.17)

v(x,t)v’T(r).

(For the sake of simplicity we write n and v instead of ñ and 5’!) With (3.17) and the derivatives formulated in T =

—~

a~,

a~ =

a~

(3.18)

the system of partial differential equations (3.2) and (3.3) becomes an ordinary one: (3.19) (3.20)

(v—T)n’+nv’=O

2n’+(o—r)v’=0, yn~ where the prime indicates the derivative with respect to

T.

The equations (3.19) and (3.20) are

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

325

nontrivially solvable, provided that (v



r)2

=

‘y n~1

(3.21)

where y n~’represents the square of the sound velocity of the system: 2

y—1

y—l

c = dpe/dne = 7 tie 7 ii = dPe/dfl• Because of the quasi-neutrality assumption tie = n, one has Pe obtains: V = T ±\/~

ne’l~ =

=

(3.22) n~.For the velocity v in (3.21) one

(3.23)

T ± C.

In (3.23) the positive sign has to be chosen to give a rarefaction wave propagating into the unperturbed plasma (negative x-direction). In the unperturbed, homogeneous plasma it should hold n = 1 and v = 0. Therefore, it follows from (3.23), that we only have to consider the solutions for r ~ The point r= dividing the regions of stationary and moving plasma, is a point of weak discontinuity. The first derivatives of n and v are discontinuous there (see [43]). Eliminating v and v’ in (3.19) by means of (3.23) and integrating, we get the density n, the velocity v, the potential 4 and the electric field E as functions of r: ~

1 ~(y+1) —

(r+~)]

2/(y—1)

(3.24)

v(r)=~i(r+\~y) ~(r)

=

y

~

(3.25)

1) ~ +

-

-

~

(3.26)

(3.27)

The solutions (3.24)—(3.27) make sense only if the sound velocity c(T) is non-negative. We express c= ~ in terms of the velocity u [see eqs. (3.24) and (3.25)]: ~

c=V~(1— ~ From c

~),

v0.

(3.28)

0 it follows, that the velocity u has to satisfy the inequality

v ~ 2V~/(y 1). —

(3.29)

If v reaches the upper limit, the density and the electric field are zero, and the potential adopts the value —y/(y 1). This result is consistent with the boundary conditions in the expansion region. —

326

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

Altogether

T

has to be restricted by the following inequalities: (3.30)

2\1~/(yl).

T~

The equations (3.24)—(3.27) differ from those of ref. [53] through the choice of the sign in eq. (3.23) only, which simply accounts for the opposite direction of the expansion. In our case at fixed t the velocity should increase with increasing x. In the limit y = 1, the set of equations (3.24)—(3.27) reduces to the well-known self-similar solutions for the isothermal expansion (Gurevich et al. [33, 43], Allen et al. [36]): n(r)=eT~

(3.31)

v(T)=T+1

(3.32)

4(r)

=

—(T

+

1)

(3.33)

with T=x/t and —1=
E=—a~4=1/t.

In the isothermal case, ‘y = 1, the electric field is space independent and does not vanish for x + so; in addition, it has a singularity at t = 0. If one initially assumes a density step, ions being located at the discontinuity x = are accelerated to infinity by the infinite electric field as soon as t> 0. The ion density decreases exponentially for t> 0 and extends to infinity. The ion velocity increases linearly in x for a given t and, hence, reaches very large values for large x (v —÷ +so for x—~+so). The rarefaction front, on the other hand, propagates with ion sound speed towards the unperturbed plasma. In the isothermal case the self-similar solution does not exhibit an ion front, and ions initially located near the density step are accelerated according to [36, 46]: —*

=

ln(t/t~)+

00

,

o,~= const

.

(3.35)

Numerical investigations of the plasma expansion problem with charge separation, where indeed an ion front occurs, claim that the time behaviour of the ion front velocity can be described by eq. (3.35) (e.g. [46]). We shall come back to this point in section 10. We terminate this subsection with some remarks concerning the validity of the self-similar theory and the corresponding solutions (3.24)—(3.27), respectively (3.31)—(3.34). In general, the existence of self-similar variables implies the lack of characteristic lengths and times [53, 54]. In the plasma a corresponding characteristic length is the Debye length. The quasi-neutrality assumption, Tie = n, being imposed at the beginning of this subsection, takes care of the fact that the Debye length looses its importance as a characteristic length. By this way, it was possible to investigate the self-similar behaviour of the plasma expansion problem. On the other hand, the group theoretical ansatz for the generation of self-similar variables can be considered as a mathematical tool to reduce partial differential equations to ordinary ones without taking care of characteristic lengths and times. This has been done by Shen and Lonngren [44] (see also [61])for the isothermal plasma expansion with charge

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

327

separation where the Debye length comes into play. In contrast to the quasi-neutral case the resulting system of equations can no longer be solved analytically; solutions are obtained only under very restrictive assumptions [48]. In any case, the existence of self-similar variables has to be checked by means of the group theoretical formalism. In cases where physically relevant self-similar solutions are found, they usually describe the asymptotic behaviour of an extended problem (see [54,61]). In the next subsection and in section 10 it will be shown that this statement holds for the plasma expansion problem, too. Next we turn to the theory of simple waves and of Riemann invariants, and discuss the properties of solutions which contain the self-similar behaviour of the plasma flow as a special case. We no longer insist on the quasi-neutrality assumption which was necessary to get the results in this subsection. 3.2. Simple waves and Riemann invariants In order to formulate the plasma expansion problem in analogy to the description of one-dimensional propagating waves of arbitrary amplitude in ordinary hydrodynamics (compare, e.g. [53]), we require that the velocity, the density and the potential are functions of each other. It is, furthermore, assumed that all quantities are piecewise monotone. Under these assumptions the equations (2.22) and (2.23)

a,n + a~(nv)= 0, a1v + o a5v

=

transform to

an

d(nv) an

at

—+

av

dn

/

—=0 ax

(3.36)

d4\au

(3.37)

We do not have to consider that the potential In view of

an/at

ax

an/ax

at

and

t~

follows from Poisson’s equation.

au/at

ax

au/ax

at

it follows from (3.36) and (3.37):

ax/atI~ d(nu) /dn =

=

ax/at~0 u + d4/du. =

u

+

n dvldn

(3.38) (3.39)

The density n uniquely determines u. It does, therefore, not matter whether the derivative is taken at constant n or u, so that

ax/atI~ ax/at~~ .

(3.40)

328

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

From the equations (3.38) and (3.39) we get n dv/dn = d4/dv

=

(d4ildn) (dn/dv)

(3.41)

and

dv/dn = ±(1/n)\/~~7’.

(3.42)

Analogous to ordinary hydrodynamics [53]we identify the expression \/~d4.ildn in (3.42) with a sound velocity c, so that cas\/~d4/dn.

(3.43)

In the following we shall use the term “pseudo-sound velocity” for c in (3.43) to distinguish it to the sound velocity in an ideal neutral gas. The relation to this sound velocity follows from (3.43) as a special case by assuming quasi-neutrality, ne = n, and by differentiating eq. (2.25) which is resolved for 4 with respect to n (see also eq. (3.22) and ref. [53]): c~:=~j’~

(3.44)

2~ ~(~‘)‘

In the case of isothermality, y = 1, CON becomes independent of the density, e.g. unperturbed, homogeneous plasma, where the ion density equals unity, it holds: CQN(n—1)V~

CON

1. In the

(3.45)

.

Assuming that the potential ~ is dependent on x only through the ion density n, one obtains from (3.43) the actual local sound velocity of the system: 2

C (x,

t) n a,b/ax an/ax =

=



nE an/ax

(3.46)

Equivalent to the above mentioned formulation of the plasma expansion is the replacement of the electrostatic field E by a pseudo-pressure J3’ in the ion momentum equation (2.23)

E=—1a~p=—a~~.

(3.47)

Using the well-known definition for the sound velocity c2 C2=j~=

dn

ak/ax an/ax

as dp/dn

(see [53]), we obtain from

(348) .

and from (3.47) again the equation (3.46) (see also [56]). Returning to eq. (3.42), its integration with respect to n orp o=±J~dn=±J n

~ n(p)c(p) -

,

(3.49)

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

329

yields a relation between the velocity and the ion density or the pseudo-pressure, respectively. From (3.47) and eq. (3.49) which is differentiated with respect to it follows: ~,

(3.50)

By means of (3.50) we substitute d4i/du in eq. (3.39) and integrate over t: x = [v

± c(v)]

t + x0(u).

(3.51)

The function x0(v) is an arbitrary function of the velocity being determined by the initial conditions; c(u) is given by eq. (3.49). The equations (3.49) and (3.51) represent the general solution of the problem and correspond formally to the solutions for propagating simple waves of arbitrary amplitude in an ideal gas (Riemann (1860), see also [53]). Both signs in (3.51) describe waves, moving relatively to the medium in the negative or in the positive x-direction. In contrast to the gas dynamic case, the u-dependence of the sound velocity can no longer be described analytically in the case of charge separation (see eq. (3.49) and (3.50), respectively). In the gas dynamic (quasi-neutral) case, we obtain for the u-dependence of c (see [53]): c(v)~(1±~v).

(3.52)

In the contrary to section 3.1, where we have investigated the self-similar behaviour of a quasi-neutral plasma, the xt-dependence of the velocity v and, therefore, of the other quantities, too, is only given implicitly. The self-similar solution is obtained as a special case of a simple wave, in which x0(u) in (3.51) vanishes (see (3.51) in comparison to (3.23) and (3.25)). The velocity of the points of the wave profile described by (3.49) and (3.51) u=u±c

(3.53)

has to be understood as the superposition of the propagation of a disturbance with sound velocity relative to the medium and of the motion of the medium itself with the velocity v. Because of the different local velocities, the profile can change its shape in the course of time. If there is a compressional wave, i.e., a certain region in which the velocity v, the density n and the pressure p are decreasing in the direction of propagation, the wave profile will progressively steepen. Finally, the wave profile can bend out of shape in such a way that the dependent flow variables (e.g. n, u) are no longer defined uniquely in space for fixed t. Three different values of n and v, respectively, are belonging to the same point x. Such a process is called “wave breaking” and opposes the foundations of the hydrodynamic description. In fact, there are strong discontinuities and shock waves, respectively, in the flow which are characterized by singularities in the derivatives and by jumps in the flow variables. Except these points and surfaces of discontinuity, the space dependence of the flow profile is unique and can be described further on by the equations for an ideal gas (see [53]). In the case where surfaces of discontinuity are formed in a flow, the wave is reflected from these surfaces. Therefore, the wave is no longer propagating in one direction, and the unique dependence

330

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

between n, u and j3’ in the equations (3.49) and (3.51) is lost. This means that generally the solution in the form of a simple wave is only valid for a finite time interval. In the special case of a rarefaction wave, however, in which the density increases monotonically in the direction of wave propagation, the solution is valid for all times. Such a wave is obtained, e.g., when a piston is drawn out of an infinitely long tube filled with gas [53]. The occurrence of the discontinuity is connected with a singularity in the first derivative of the velocity. In the case where the system does not contain any additional internal boundaries, the second derivative also becomes singular. Therefore, it is possible to determine the instant and the location at which the wave ceases to exist. For the inverse function x(v, t) it must hold [53]:

a2x/av2j,

ax/av~,= 0,

=

0.

(3.54)

For the gas dynamic (quasi-neutral) case, where c(v) is given by (3.52), one obtains by differentiation of (3.51) with respect to o from the conditions (3.54): ~

2 y+1 av



av2

~

(3.55)

The steepening time t

5 depends on the polytropic exponent y, and on the slope of the function x0(v) obtained from the general solution of (3.51) where ax0/av <0 represents a compressional wave. Each compressional wave steepens in the course of time and develops due to (3.55) at the time t~a vertical tangent in the velocity profile v(x). In ordinary hydrodynamics the presence of such discontinuities in an otherwise ideal flow produces an energy dissipation which causes a strong damping of the wave. The strength of the dissipation can be determined by applying the conservation laws of mass, energy and momentum to both sides of the discontinuity yielding the Rankine—Hugoniot equations [53]. Furthermore, the appearance of dissipative effects causes an entropy increase and, therefore, irreversibility. In the plasma, however, dispersive effects arising from charge separation may counteract the progressive steepening without affecting the reversibility. Montgomery [34], who investigated nonlinear ion acoustic waves, suggested that the steepening of the wave can be balanced by dispersion. Charge separation effects come into play when the scale length of the electrostatic potential is comparable to the Debye length. This suggestion is certainly true for small amplitude perturbations and weak inhomogeneity, respectively (see Sagdeev [55]). However, if there is an extremely inhomogeneous and nonlinear expansion of a plasma into vacuum, which is typical for the initial stage, unlimited steepening and the resulting wave breaking will occur despite dispersion. This will be shown in sections 5 and 6. An indication for this statement is the general formulation of the expansion problem with charge separation as given by (3.49) and (3.51). In contrast to the quasi-neutral case [see (3.55)], the steepening time t5 cannot be calculated in a simple manner, since c(v) is not given analytically. A further relation between the plasma expansion problem and ordinary hydrodynamics results from the construction of Riemann invariants (see [531).For this purpose we thewhere following 2/n)use a~n, c2 is trick givenand by replace For the the electrostatic field ion in the ion momentum by obtain: —(c (3.46). hydrodynamic equations (2.22) andequation (2.23) we a,n+oa~n+na~o=0

(3.56)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

a1u+ua~u+Ea~n=o.

331

(3.57)

By multiplying eq. (3.56) with ±cln,and by adding it to (3.57), one obtains: (3.58)

[a1+(u±c)a5]R+=0 with

f

R ±:= v ± dn’ n’c(n’) = J

~ Jf c(p’)dp’n(p’)

(3.59)

±

R~and R are the Riemann invariants. The differential operators acting on R~and R turn out to be the derivatives in the direction of the characteristics C~and C in the xt-plain, given by (dx/dt)~= u

+

c,

(dxldt)_

=

u



c

(3.60)

(see also (3.53)). On each characteristic C.. and C, respectively, the corresponding Riemann invariant remains constant. This kind of formulation represents an alternative to the numerical treatment of the expansion problem by finite difference methods. By means of the Riemann invariants the hydrodynamic ion equations (3.56) and (3.57) can be formulated as a Cauchy initial value problem, and can be solved by the method of characteristics. In the quasi-neutral, isothermal case (y = 1), where R.,. is given by R~=u±lnn

(3.61)

this was done by Felber and Decoste [9]. The comparison with eq. (3.49) shows the correspondence between Riemann invariants (3.59) and simple waves. With R~asO we obtain from (3.59) —

f

dn’ c(n’) n’

from which follows due to (3.46) a~u=~a~n=±~E. The replacement of E in the ion momentum equation (2.23) by

[a1+(u~c)a~]u=0.

(3.62) ±c a~vyields

finally: (3.63)

With (3.63) we thus have succeeded in formulating the ion momentum equation (a1 + u a~)u= E in the form of a simple wave equation. The pseudo-sound velocity c(x, t) defined by (3.46) is introduced by the analysis of simple waves of arbitrary amplitude [see (3.49) and (3.51)] or by the construction of

332

Cli, Sack and H. Scharnel, Plasma expansion into vacuum — A hydrodynamic approach

Riemann invariants [see (3.58) and (3.59)]. The pseudo-sound velocity can be obtained as a function of x and t by evaluating the numerical data for the ion density n and the electrostatic field E, the latter following from Poisson’s equation. It can also be used to characterize the state of charge separation. We shall come back to this point in section 5. Furthermore, c(x, t) from (3.46) allows a conventional stability analysis of the numerical difference scheme (see section 8 and appendix D). In the next section we discuss the foundations of the numerical methods for solving the hydrodynamic ion equations (2.22), (2.23) and Poisson’s equation (2.24). Of special importance will be the derivation of the boundary conditions adopted to the finite space interval.

4. Numerical methods and analysis 4.1. Lagrangian equations of ion motion In many hydrodynamic problems, especially in one-dimensional systems, it is useful to transform Euler’s equations of motion to a coordinate system moving with the local flow velocity. In this so-called Lagrangian picture the equations of motion get a very simple form, which can be handled more easily, analytically and numerically. Often the analytical investigation of hydrodynamic problems is not possible until Lagrangian coordinates are introduced [62, 63, 64]. In section 6 we shall discuss the consequences of such a formulation in the concrete case of the plasma expansion. For the derivation of the Lagrangian equations from the Eulerian picture we assign a certain label ~ to each fluid element; ~ is the Lagrangian space-coordinate and represents the position of a fluid element at the time t = 0. The relation between the space-coordinate x and the local flow velocity u is given by the equation of motion dx(t) /dt = u(x(t), t)

.

(4.1)

The integration of eq. (4.1) yields implicitly the coordinate x of the fluid element for t >0:

x(t) = ~ +

v(x(T), T) dT

(4.2)

with v(x(r),r)=v(~+Jo(x(r’),r’)dr’,r)as5’(~,r).

(4.3)

From the equations (4.2) and (4.3) we obtain the transformation of the Eulerian variables (x, t) to the Lagrangian ones (~, T) [42, 50, 56, 62]: x=x(~,r)=~+Jff(~,r’)dr’

(4.4)

t=t(~,r)=r.

(4.5)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

333

According to (4.4) and (4.5), respectively, the partial derivatives a1 and a~transform as follows:

a

a

a

d (4.6)

+

= (i

a5’(~,r’)dr’)

(4.7)

~.

Equation (4.6) defines the convective time derivative. The transformation (4.7) is uniquely determined only if

J

a5’(~’)dr’

-1.

(4.8)

The equality sign means that in the xt-space two particle trajectories cross each other. In this case the dependent variables are no longer unique, and the plasma cannot be described by ideal hydrodynamics (see sections 5 and 6). Inserting the equations (4.6) and (4.7) into the ion continuity equation (2.22) and the ion momentum equation (2.23), we obtain together with eq. (4.1) the description of the ion motion in the Lagrangian formulation [42,50,56]: x(~,t)

=

n(x(~,t), t) v(x(

~,

u(x(~,t), t)

(4.9)

n~(~)

(4.10)

=

~

t), t) = E(x(

~,

t), t)

(4.11)

with the initial conditions given by n(x(~,0),0) = n0(~)

and

x(~,0)=

For solving the continuum problem numerically, the equations (4.9)—(4.11) are discretized in space and time [65,66]. An essential advantage of the Lagrangian formulation is the absence of the nonlinear convective term u a~u,which causes most trouble in the Eulerian picture. The equations of motion can easily be discretized via an explicit Lagrangian difference scheme. This procedure excludes numerical dissipation and dispersion which usually lead to an undesirable strong smoothing of the solution. Such a falsifying influence must be expected if the Eulerian equations of motion are discretized by an explicit Lax scheme [65,66]. The Lagrangian method also guarantees that the difference mesh is coupled to the local flow velocity, so that strong inhomogeneities in the velocity and in the density profile are resolved’ in an optimum manner. Further details of the discretization procedure are to be found in appendix C. The explicit difference scheme turns out to be a special case of the more general formulation. For solving the equations (4.9)—(4.11) the knowledge of the potential ~ and of the electric field E, respectively, is required. In the next subsection we present two solution methods for the nonlinear Poisson’s equation.

334 4.2.

Ch. Sack and H. Schamel, Plasma expansion into vacuum



A hydrodynamic approach

Two iterative methods for the nonlinear Poisson’s equation

After replacing the electron density by the electrostatic potential, Poisson’s equation becomes a nonlinear differential equation [except y = 2, see eqs. (2.24) and (2.25)1. Its numerical solution causes some problems with respect to the numerical procedure and to the formulation of the boundary conditions. Because of the nonlinear term [1 + ((y 1)Iy)~]’~’~’~, y = 1 and ed’, y = 1, respectively, the electrostatic potential is not available by a direct spatial discretization of Poisson’s equation (2.24), even if the ion density is known. To get rid of this difficulty, it is necessary to develop suitable iteration procedures. In the first iteration procedure used by Mason [38], the potential 4 is replaced by the electric field E = —4’. Differentiating Poisson’s equation (2.24) with respect to x, and expressing the electron density tie by n E’ yields a differential equation of second order for the electric field: —



(4.12)

E”—~(n—E’)2~E—n’=0.

In the isothermal case, y = 1, eq. (4.12) simplifies to that equation solved by Mason [38].Insofar (4.12) represents an extension of the model to more general polytropic equations of state. In implementing the iteration procedure, we assume that the electric field E is known for the iteration step p at each grid point i. In order to obtain E for the next iteration step p + 1, we insert E of the pth iteration step at the grid points i + 1 and i 1, and solve with respect to E,. The advantage of this procedure is that we get E at each grid point i by solving a simple algebraic equation. The details of this iteration scheme, especially the discretization of (4.12) on a non-equidistant grid, have been described in ref. [56] for the isothermal case, y = 1. The second iterative procedure, which yields the potential as the solution, used the Newton method and was developed among others by Dewar et al. [67] and Forslund et al. [68]. It is based on a Taylor expansion of the nonlinear electron density under the condition that 4 can be described by 4 = + ~. Linearizing Poisson’s equation (2.24) with respect to gives —

~,

dn

(4.13)

t~~——4~+n(4~)—n

c~’~ ‘I’D

with —

1

lI(y-l)

Tie(~o)= (i+ ~

y~1

(4.14)

y1

and 1

dn d~

=

0

-

~

(i + 7



1

1

(2—y)/(y—l) =

-

[ne(~o)]2~,

~

~ 1~

(4.15)

y1

A version which can be handled numerically more easily, is obtained by subtracting the term

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

335

(dneldq5o) q50 from both sides of (4.13): dn

dn (4.16)

In contrast to (4.13), eq. (4.16) contains only one spatial derivative of second order. In addition, eq. (4.16) represents the more suitable form to implement the boundary conditions which are formulated in q~instead in 4’,. The spatial discretization of eq. (4.16) yields a tridiagonal coefficient matrix for the potential 4’ at the grid point i, the first and last row of which are filled up by the boundary conditions (see section 4.3). A detailed description of the discretization procedure on the non-equidistant grid will be given in appendix B. Prescribing the initial value 4~ and the boundary conditions, we obtain a preliminary solution for 4’, which is inserted instead of 4~in the next iteration step. In the case of convergence, i.e., 4’ 4.~—>0, one obtains from eq. (4.16) the original Poisson’s equation (2.24). For t = 0, the quasi-neutral solution is used, given by —

~‘

4’0(x,t=0)=

y—1 lnn,

(n~—1),

y~1 .

(4.17)

y1

It is found, however, that the iteration procedure is rather insensitive to the choice of the initial function. Test calculations with 4’0(x, t = 0) as 0 have given identical solutions compared to those with (4.17) as initial function. However, more iterations have been necessary to gain convergence. Therefore, eq. (4.1) has been used as initial function in all of the calculations within the second iteration procedure. For t > 0, 4~is identified with 4’ from the preceding time step as the initial function for the iteration. An essential advantage of the second iterative procedure is the fact, that besides 4” further spatial derivatives as, e.g., n’in eq. (4.12) are not involved. Thus, the convergence is not deteriorated by an inaccurate calculation of possibly steep gradients and of small scale structures. In section 5 we shall see that this property of the iteration procedure plays an important role for the resolution of nonlinear structures occurring in the expansion problem. Moreover, the second iteration procedure converges considerably faster than the first one. After ten iterations, e.g., the maximum relative error (see appendix B) amounts to about 108. Using the first iteration procedure, many more iterations are needed to reach this degree of accuracy [56], especially in the region of low density. Because of these properties we prefer the second iteration procedure, henceforth called “Poisson solver”. We note that the second iteration procedure can be applied to any equation of state for the electrons [68]. In the present case, only those equations of state are considered which imply an electron front at infinity. In section 2.2 we have shown that this is the case for 1 ~ y 2. For y >2 (and for an equation of state describing trapped electrons, respectively), dfle/d4’0 in eq. (4.16) becomes singular when the electron density ne goes to zero [see eq. (4.15)]. In this case, the Newton procedure cannot be used appropriately. A way out of this dilemma is offered by the fitting procedure of Crow et al. [42],where the ions are given as a step function. In this case a numerical solution of Poisson’s equation in the region of the pure electron cloud is not necessary since there exists already an exact analytical solution (see section 2.2). Therefore, Poisson’s equation has to be solved numerically only in the region where

336

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

the ion density is non-vanishing, for which the second iteration procedure can be used again. The resulting solution is then matched to the analytical solution from the pure electron region. In this and in the last subsection we have discussed the basic properties of the procedures for solving Poisson’s equation and the hydrodynamic ion equations. A complete solution of the equations (4.9)—(4.11) and (4.16), however, is not possible until the general boundary conditions, eqs. (2.27), are adapted to the numerical solution method. In the next subsection we shall derive such adjusted boundary conditions. 4.3. Boundary conditions It is clear that numerical solutions can only be obtained in a finite integration interval. In order to solve systems with boundaries at infinity, one has to transform them to a finite interval. In the present case, the length of the integration interval is 2L, L ~ x c + L, at the beginning. The boundary conditions for this finite interval at x = ±L have to be formulated such that the boundary conditions of the continuous problem are satisfied automatically, when x ±so~This concept leads to so-called open boundary conditions, which are represented by differential equations. Thereby it is not the function itself which is prescribed at the boundaries, but its asymptotic behaviour in differential form. Firstly, we investigate the asymptotic behaviour of the system (2.22)—(2.25) formulated in Eulerian form for the undisturbed, homogeneous plasma, x—~—so and x = —L, respectively. It is assumed that the density, the velocity and the potential vary only weakly in this region so that linearization of the system with respect to small perturbations is justified. We write —

—>

n=n0+5’,

4’=4’~+4,

v=u0+5’,

ne=neo+ne,

(4.18)

where n0 = neo = 1, 00 = 0, and 4’~= 0 are given by their values at x = —so [see boundary conditions in section 2.1, eqs. (2.27)]; the perturbed quantities are marked by tildes. With (4.18) we derive from the equations (2.22)—(2.25) the linearized system of equations

a1i~+a~5’=0

(4.19)

a,5’= —a~4’

(4.20)

1-cb—n.

(4.21)

a~=

Making a Fourier ansatz f=~f~exp[i(kx—wt)] with

(4.22)

Jss (,~,5’, ~)

and

fk as (ilk’

Uk, ~k)’

and inserting it into (4.19)—(4.21), we obtain the well-known dispersion relation for ion acoustic waves: 2 (0 k

=

yk 1+

2 7k

resp.

~k

±_÷ ——

~f~k

\/1+

2 7k

wkE’~.

(4.23)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

As seen from the definition, w~and

Wk

337

differ in the sign in front of k, i.e., w(w~) represents waves

propagating to the right (left), (4.24) Demanding that no waves should come in from x = the electrostatic potential

4’ =

~

k>0

~ exp{i(kx + wkt)} + k<0 ~

with w~>0 for k >0, and

Wk

—so

at x =

—L,

we obtain for the perturbation

~k exp{i(kx + w~’t)}

(4.25)

<0 for k <0. From the reality property of

4’,

it follows (4.26)

= 4’k~

~

where use is made of Wi~ik= Equation (4.25) then becomes

=

~

~k

4’ of

exp[i(kx

+ ~

= ~

~k2

and asterisk means complex conjugation.

~)]

(4.27)

which in the long wavelength limit (k ~ 1) reduces to ~—~4’kexp[i(kx+v5~

(4.28)

kt)].

Because of k ~ 1 the dispersive term being proportional to k3 has vanished. (It is well-known that this term represents dispersion in the Korteweg—de Vries equation [62] which is derived for weakly nonlinear waves.) From (4.28) we immediately get by differentiation the corresponding differential equation we are looking for:

a,çb where



4’ is

\/~

a~cb= 0,

replaced by

4’ [see

(4.29) (4.18)]. With E =

—4” we also

get

a

1E—\/3a~E=0,

(4.30)

and the asymptotic behaviour of the velocity in the unperturbed, homogeneous plasma is described, of course, by a differential equation of the same kind: a1u —\/~ a~v=0.

(4.31)

The space-time discretized version of (4.31) is used to determine u at the left boundary x = L in the numerical solution (see appendix C). Thus, it is not the value of v that is prescribed but its variation, in differential form. One has not to distinguish between the Eulerian and the Lagrangian picture, since the term u a~vis lost after linearization. —

338

Cli. Sack and H. Scliamel, Plasma expansion into vacuum — A hydrodynamic approach

In principle it is possible to use (4.29) as a boundary condition for Poisson’s equation at x = L. Since in Poisson’s equation time appears only as a dummy variable, we have chosen another boundary condition containing only spatial derivatives. This so-called “gradient boundary condition” has already been derived by Crow et a!. [42] for the case of the isothermal plasma expansion. We demand again that the density and the potential at the left boundary, x = L, vary weakly. Furthermore, we assume that the variation of the ion density can be described by a power series in 4’ with time dependent coefficients a,,, v >0, ii E ~J: —



n=1+E,

(4.32)

e=~a,,4’~.

Poisson’s equation linearized with respect to

~



1

4’)

l/(y—l)

4’

(2.24) then reads:

1 1 —n1+—4’+...—(1+ai4’+...)~(——a1)4’.

After multiplication of (4.33) with x—> —so we obtain:

4”

and integration using the boundary conditions

(4.33)

4’ = 4” = 0 for (4.34)

Since the right-hand side of (4.34) is constant in x, we can write (4.35)

= 4”4’~-L+~~,

where z~xdenotes the spatial step. In a similar manner we obtain the gradient boundary condition for the E”-equation (4.12) (see also ref. [56]): E’/ELL

(4.36)

= E’/ELL÷~X.

From the spatial discretization of (4.35) and of (4.36) we get for 4’ and E the left-hand boundary condition. It depends on the choice of the iteration procedures described in section 4.2. Note, that in the contrary to (4.29) and (4.30) the relations (4.35) and (4.36) are nonlinear. Consistent with the corresponding iteration procedures, the boundary conditions (4.35) and (4.36), respectively, are iterated as well. Hence, the behaviour of the plasma quantities at the left-hand boundary is completely determined by (4.31), (4.35) and (4.36), respectively. In the vacuum region, x—> + so and x = +L, respectively, the boundary conditions must be determined separately. Again it is our intention to derive differential equations of the first order, describing the asymptotic behaviour of the solution for x + so correctly. A differential equation of this kind for the potential has already been derived in section 2.2. This equation resulted from the integration of Poisson’s equation by neglecting the ion density: 2 1 4’) y/(y—l) =~=n~ E ~=(i+~ (4.37) —>

4’F2



Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

with

Tie

from (2.25). For 2. = —“/~ n,,~’

4”

339

it follows from (4.37): (4.38)

A further integration of (4.38) to determine the boundary conditions is not meaningful because two integration constants, xf and 4’~,are introduced, which are not available until the solution of the entire set of equations is known (see section 2.2). According to the preparation of the second iteration procedure, eqs. (4.13)—(4.16), the nonlinear term n~2in (4.38) is linearized by a Taylor expansion: n~2asF(4’)~F(4’ 0)+~4’i;

4’=4’~+4’~,

4’1HI4’0~

(4.39)

with 2 and dF/d4’ [n,,(4’0)]~” 0= 2from (4.39) into (4.38) and add on both sides the term (1/V’~)[ne(4’o)]1’~”24’o: We insert n~ F(4’0)

4” +

~

=

[(4’~)]~4’

=

~

{[n~(4’~)]~ 1[ -

(4’)]1~2

4’o}.

(4.40)

The spatial discretization of (4.40) determines the boundary condition for Poisson’s equation at x = +L. Consistent with the difference scheme for (4.16), eqs. (4.35) and (4.40) fill up the first and the last row of the tridiagonal matrix, the solution of which yields the potential 4’ at the grid point i [see appendix B]. Using the first iteration procedure, we obtain the corresponding boundary condition for x = + L from the integration of the E”-equation (4.12) with n = n’ = 0:

E’ + (E2I2)1~= 0. The velocity v at the right-hand boundary, x = [see appendix C]: du/dt=E.

(4.41) +L,

is determined by the spatial average of the equation (4.42)

The boundary conditions at x = ±L, are, therefore, given completely by eqs. (4.31), (4.35) and (4.40), (4.42), respectively. In the next section we formulate the conservation laws derived from the set of equations (2.22)— (2.25), which turn out to be important tools for testing the reliability and the accuracy of numerical solutions. 4.4. Conservation laws Multiplying the continuity equation (2.22) with u212, and the ion momentum equation (2.23) with no, and adding the resulting equations, one gets in the Eulerian description

a,(nu212) + a~(nu3I2)= nuE.

(4.43)

340

Cli. Sack and H. Schamel. Plasma expansion into vacuum — A hydrodynamic approach

In contrast to an opposite statement made by Mora et a!. [47]and by Trulsen [69],it will be proven that eq. (4.43) can be formulated in conservative form. Eliminating no through tieVe — a~Eby Ampere’s law, and using the electron continuity equation, one gets

a1(nv2/2 + E2/2) + a~(nv3/2+

4~atne

4’neoe) =

(4.44)

.

It is worth noticing that this equation holds true even if the electron inertia is included. In the inertialess limit the right-hand side of (4.44) can be expressed as a partial time derivative

—4’ a,n~= —atFe(4’),

(4.45)

where Fe(4’) =

ne(4’)

noting that p~(4’)=

4’



(4.46)

Pe(4’),

tie(4’)

and ne(4’) is given by (2.25). We therefore, arrive at the following energy

conservation law

a,[nv2/2

+

E212 + Fe(4’)]

+

a~[nv3/2 + (neve)4’]

=

0,

(4.47)

which is in conservative form. The quantity ~‘,, defined by

=

J

dx [nv2/2+ E2/2 +

(4.48)

Fe(4’)],

does not change in time. Another conservative formulation of the energy law can be found. First, we notice that Fe expressed in terms of Pc becomes in the limit of inertialess electrons

Fe(Pe)=

Pc 7 1/y y1 y—1 Pc p~(lnp~~l), —



y1 7=1

,

(4.49)

where use is made of p~= n~and of eq. (2.25). With this form of Fe, eq. (4.47) transforms to +

ç

+ y~1)

with 4’~= 4’(x = so, t) The quantity ~‘2

~2

=

+

—y/(y

Jdx[no2/2+E2/2+

a~[~+ —

Tie Ve

(4’~&)] =0

(4.50)

1), y ~ 1 (see (2.27)).

~],

(4.51)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

341

is constant in time as well. It represents the total energy and reflects the balance between ion kinetic energy, electrostatic field energy and internal energy of the electrons. On the other hand, taking into account electron inertia effects, it follows from the electron momentum equation (see eq. (2.16)) = fle[1 +

where 6

=

6E’ (arve +

Oe

axve)],

(4.52)

me/mi. Equation (4.44) then becomes

a[ç +

+

Fe] +

a5[~ +

1 fleVe4’]

+

6[TieE

a4’ (a

1v~+ o~axue)]

=

0.

(4.53)

In the inertialess limit, 6—*0, we recover the old expression (4.47). A straightforward calculation shows that (4.53) is equivalent to the following conservative formulation:

a,[ç+6

~.~+ç+ ‘~1]+a1[~+o

~-~+

~‘1 VePe]reO,

(4.54)

which is the well-known expression for the energy conservation using a polytropic electron equation of state. This equation reduces to (4.50) in the limit 6—*0 using Pc n~and 7

7

y—l VePe= )‘l

/ neuel~1+

y—1\ y

Similarly, the momentum can be expressed in conservative form, 2 + Pc E2/2] = 0, a,(nv) + a~[nv where E = —~~4’ and Poisson’s equation has been used to replace —n a~4’by a~[E2/2 PeL —

(4.55)



The total momentum p= fnudx

(4.56)

increases linearly in t: P(t) = t,

(4.57)

which can be seen by calculating P(t) and using the boundary condition p~—~ 1 for x—> —so, noting that the other quantities vanish at the boundaries. It is the electron pressure at minus infinity which pushes the ions and gives rise to a monotonic increase of the total momentum. The conservation laws (4.47) and (4.50), respectively, and (4.55) are in the appropriate form for implementation in the numerical code.

342

Cli. Sack and H. Scliamel, Plasma expansion into vacuum — A hydrodynamic approach

4.5. Variational principle The energy conservation law (4.47) suggests the following Hamiltonian H=

J

dx [no2/2 + (a~4’)2/2+ Fe(4’)]

(4.58)

for the formulation of a Hamiltonian principle. Introducing the potential representation of the curl-free velocity field (Clebsch [70], Elsässer and Schamel [71]) o

—asp

=

(4.59)

,

we get by variation with respect to ~, n and 4’ under the constraint that the 4’ variation is dependent on bn and satisfies Poisson’s equation, the following Hamiltonian set of equations

a,cp=~H/~n= v2/2+ 4’, a

1n

=

(4.60a)

—~HIF~p = —a~(nv).

(4.60b)

The first one is the well-known Bernoulli equation, that is the integrated form of the momentum equation, and the second one is the continuity equation. We note that the Hamiltonian density ~1C in (4.58) or the equivalent formula in which Fe(4’) is replaced by (4.49) can be written up to a constant as

2/2 + ne W(ne)],

=

[n02/2 + E

(4.61)

where

—J

W(n~):=

(Pc



1) d(1/fle),

(4.62)

is the compressional energy of the electrons per unit fluid mass (Elsàsser and Schamel [71]). The Lagrangian density from which the Hamiltonian density can be derived reads ~t= n[~

-

~(a~~)2] ~(a~4’)2 F~(çb), -

-

(4.63)

from which immediately follows that ~ and n are canonically conjugated variables pas6L/&~b=n,

(4.64)

being already implied by (4.60). The Euler—Lagrangian equations minimizing the functional j’ dt J dx ~ are again given by the set of equations (4.60a, b) and (2.22), (2.23), respectively, where the ~4’variation is, as before, related to ~n under the constraint of Poisson’s equation.



Ch. Sack and H. Schamel, Plasma expansion into vacuum A hydrodynamic approach

343

Hence, we see that our basic equations (2.22)—(2.25) can be completely derived from a variational principle. In the next section we present the numerical solutions obtained by the dissipationless numerical scheme described in this section. Our emphasis is to clarify the role of charge separation, and of charge neutrality, respectively, and of the electron thermodynamics on the dynamical evolution of the system. In view of the intimate relationship to ordinary hydrodynamics as discussed in section 3, it will not be surprising that almost all solutions exhibit discontinuities. We shall describe these discontinuities and shall offer ansätze to understand the singular behaviour of the plasma flow.

5. Numerical results for the dissipationless plasma expansion Figure 1 shows the space-time dependence of the ion density in the isothermal case, y = 1, taking into account charge separation. The global behaviour of the ion density for three different times is shown in fig. la. The framed part in fig. la is drawn on an enlarged scale in fig. lb for five successive time steps. Starting with a diffuse ion density profile (1 = 4 in (2.28)), a sharp ion front quickly develops in front of which a pure electron cloud exists. Passing a plateau-like state, a spike is formed at the front growing explosively. Finally, at t = t~, 17.94, the numerical solution breaks down without any indication of a numerical instability (e.g., oscillations from grid point to grid point). In section 8 we shall describe such a case of numerical instability. Neither the increase of the number of grid points nor the introduction of an ion pressure and of various artificial viscosities could avoid the singular behaviour of the plasma flow within an explicit discretization scheme (see also ref. [56]). Figure 2 gives a first information about the origin of the singularity. In fig. 2a the ion density n, the electric field E and the ion velocity v of the initial state (v as 0 for t = 0) are compared with the n(xt) tO

ri(x,t) 0.15

0.1

~

::

a

\\

t~13 ,‘~

\\z-Z~ \‘\—‘ r._

0

~

~~~•-~•ó

-,,~— Z~.A7

~

Fig. 1. Ion density n as a function of space and time for y = 1 in the case of charge separation; (a) gIooal behaviour of the densityfor three different time steps; (b) expansion front on an enlarged scale, —20~ x < 10.

344

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach 1.0

0.75

v.25

_‘_\

t0

02

n

E

Q15

:i

0.5

0.1 025

‘~

Q05

a O ..,.,. -200

: ~•1~- ,‘‘..

.

-100

to

0

\

/1

I

t=17

0 100

Q25 3°

iL v

02

11

I\ \/i~. / t’t~.

\I ~

E 20 0.15

~~

0.5

\/

~\

0.1

S

025

0 -200

DOS

.,,~,....,..,.

-100

0

100

0

Fig. 2. Ion density n (solid line), electric field E (dotted line), and ion velocity u (broken line) as functions of x for y two different times: t = 0 (a) and t = 17(b) for u(x, 0) 0.

=

1 and charge separation at

corresponding quantities at a later time just before the break-down. At the beginning the electric field is hump-like. During the elapse of time the velocity profile exhibits this structure, too. In the early stage of the evolution, the velocity maximum stays a little bit behind the ion front causing the narrowing of the mesh point distances ~ as x~.÷1 x1. in the region close to the ion front. Due to the conservation of —

mass given by ni.~1 = const. in the Lagrangian picture (see eq. (4.10)), the density increases in this region. At the same time the velocity and the electric field strongly steepen (see fig. 2b). The main part of the plasma in the rear behind the front approaches the state of quasi-neutrality and of self-similarity; this is expressed in fig. 2b by the electric field forming a plateau, and by the linear velocity increase [see also eqs. (3.31)—(3.34)]. The quasi-neutral self-similar region broadens in the course of time. A more detailed picture of the solution near the ion front will be presented in fig. 7. Assuming quasi-neutrality from the outset we obtain a completely different solution, although the initial density profile is the same. Figure 3 shown n, E and u for t = 0 (fig 3a) and for t = 30 (fig. 3b). At any time the ion density decreases monotonically, whereas the velocity and the electric field increases monotonically. In contrast to the case discussed above, y = 1 and charge separation, the solution does not break down. This behaviour is due to the altered initial condition for the electric field. In the

Ch.

Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

345

to 0.3

t~O

E

02:

E

:‘

~

—200

-100

n

~

0.5

02:

100

0

V~Et

\ ~.‘

-200

-100

602

V\/!

b

0.3

.

0

100 x

Fig. 3. y = I and quasi.neutrality: t = 0 (a) and t = 30 (b).

quasi-neutral case with y

=

1, the electric field is given by

E=—a~lnn=—1a~n.

(5.1)

Taking into account the prescribed ion density profile, eq. (2.28), the evaluation of the term —(1/n) a

5n

for x —*

+

so

yields

E=1/l=0.25,

(5.2)

i.e., the electric field does not vanish in the region of low density. The velocity assumes the value v = Et (see fig. 3b). The smaller value of the plateau corresponds to the self-similar solution, eqs. (3.31)— (3.34). In contrast to the case of ‘y = 1 and of charge separation, the spatially increasing ion velocity causes the fluid particles to move apart in the low density region and, therefore, the solution does not break down. One can get a hump-like structure of E and, therefore, of v, too, if one uses a different density profile decreasing, e.g., algebraically, n x~,a >0 for x—+ +so. In this case a break-down of the solution even for y = 1 and quasi-neutrality occurs. Here, therefore, the rather unphysical ‘-~

346

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A liydrodynamic approach

behaviour of E, namely E(x—* +so) = 1/1 = const., which is due to the special choice of the initial density profile, does not allow the collapse. This is, of course, an exception. In general, one has to expect that a collisionless expanding plasma is exposed to this type of collapse, as it will be shown in the following. Maintaining quasi-neutrality, such a situation is found if one chooses a different polytropic exponent keeping the initial density profile. In fig. 4a and fig. 4b, respectively, n, u and E are plotted for t = 0 and t = 10 in the case y = 2. The initial profile for the electric field shown in fig. 4a follows from [see (2.21) and (2.28)]:

E=_2a~~)=~_sech(~-_1-~°), 1=4.

(5.3)

This initial situation is very similar to that in fig. 2a, and, therefore it is not surprising that the solution also breaks down. The collapse, however, takes place at an earlier time. As in the preceding cases, the velocity assumes the structure of the electric field from the initial state (see fig. 4a). Due to (5.3) the 1.0’ n ‘0.15 t=O

0200

-100

0

x-_::

1.0

to 1.0 E

V

0~75

08

~

08

n 0.6 ~6

Q:5~\Q2

Fig. 4. y = 2 and quasi-neutrality: t

02

=

0 (a) and t = 10 (b).

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

347

field variation is enhanced (e.g., the electric field decreases more quickly for x—* +so than in the case y = 1 and charge separation). This affects the velocity and, as a result, the collapse is speeding up. Figure 4b is of great interest because it undoubtedly proves that wavebreaking associated with ion bunching does not rely on the assumption of charge separation. It occurs in the quasi-neutral plasma as well. Since the quasi-neutral description covers at the same time also the expansion of a neutral gas, we conclude that bunching is also to be seen in the expansion into vacuum of an inviscid ordinary gas. An analogous discontinuity, as mentioned already in subsection 3.2, develops in the gas flow driven by an advancing piston. The straight characteristics C~ of a simple compression wave (sometimes also called condensation wave) originating at the piston form a cusp-shaped envelope in the t—x plane. The tip of this cusp marks the break down of the flow description through a simple wave, resulting in a discontinuity and singularity, respectively, of the flow pattern [53, 72]. A phenomenon of this type was perhaps first noticed by Stokes in 1848. To realize this breaking phenomenon numerically, a code is necessary in which dissipation is more or less absent. Dissipation usually incorporated intrinsically in the numerical codes (e.g. LAX codes) would smooth out the density peak, giving rise to the known shock-type structures. We shall come back to this point in section 9 in the description of fig. 10. For the sake of comparison, the case y = 2 and charge separation is shown in fig. 5a, t = 0, and fig. Sb, t = 12. The collapse takes place at a higher value of the density similar to the case y = 2 and 1.0

I

0.15

t~0 0.75

E :~tivE

.

01

i\I

05

005

:‘\~

a

—200

1.0

-100

0

t~12

100

1.25 0.25

___-_______~....~,,~

0.75

a

\

tO

v

,Ii~E

\j

0.2 E

0.75 0.15

0.5

0.5 ‘0.1

0.25

/

0

0.25 0.05 ‘4.—

-100

-50

0 0

Fig. 5. y = 2 and charge separation, t = 0 (a) and t = 12 (b).

348

Ch. Sack

and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

quasi-neutrality. This is due to the smaller width and the stronger decrease of the electric field at the time t = 0. Remember that the electric field falls off exponentially [see (2.39)] in the region where the ion density is negligible. In the case y = 1 and charge separation, however, the electric field falls off only algebraically [see (2.36)]. The collapse time is comparable to that for y = 2 and quasi-neutrality. Behind the ion front again a self-similar region forms. For y = 2, the electric field from the self-similar solution, however, falls off linearly in x [see (3.27)]. This decrease can be seen in fig. 5b behind the ion peak. For values of y within the interval 1 ~ y ~ 2, e.g. y = 1.2 (see ref. [57]), we obtain under the same initial and boundary conditions similar solutions. Figures 3 and 4 point out the limits of the quasi-neutrality assumption in the plasma expansion. The use of quasi-neutrality is justified only if the inhomogeneity lengths are large compared to the Debye length. This is no longer true if the density, the electric field and the velocity steepen and vary on the Debye length scale (see, e.g. fig. 4b). Through this the physical assumption of the quasi-neutral description loses its validity already before the numerical solution breaks down. The Riemann solutions for simple waves of arbitrary amplitude derived in section 3.2, express the break-down of the quasi-neutrality assumption analytically (see also [34]). Hence, it is necessary to include charge separation effects in cases where profile steepening takes place (see section 3.2). In general, the solution of the complete model including charge separation effects has to approve whether and in which way the plasma approaches the quasi-neutral state. In some cases (see fig. 2 and fig. 5), this approach can be investigated by comparison with the self-similar solutions. An alternative procedure is given by the pseudo-sound velocity c(x, t), the square of which [see (3.46)] is shown in fig. 6a for three different times as a function of space for —60 ~ x ~ +20, y = 1 and charge separation (see fig. 1 and fig. 2). The plot of c2 given by c2(x, t) = —nEIa,~n

(5.4)

has been calculated from the numerical data for the ion density n and the electric field E. Figure 6b shows n and E for t = 9, and t = 15 in the interval —40 ~ x 0. In fig. 6a we can clearly distinguish two regions. The left-hand part in fig. 6 represents the quasi-neutral region which expands in the course of time and which in the isothermal case, y = 1, is described by c2 = 1 [see (3.44) and (5.4), respectively with E from (5.1)]. As it is expected, there are significant deviations from quasi-neutrality within and ahead of the ion front region which can be seen on the right-hand side of fig. 6a (see for comparison fig. 6b). The strong increase of c2 is caused by the formation of a plateau in the ion density. For a~n—+0 and E >0, E 0, c2 from (5.4) becomes singular. When the ion hump appears, there exists a region with c2(x, t) < 0 behind the ion front (see also fig. 1 and fig. 2). In this case the pseudo-sound velocity c(x, t) is imaginary and, therefore, cannot be interpreted anymore as a velocity. Consequently the transformation of the plasma expansion problem to a hydrodynamic one (see section 3.2) using the concept of pseudo-sound velocity is not applicable globally. In particular the construction of Riemann invariants, eqs. (3.62) and (3.63), becomes meaningless. In those regions, however, where c2(x, t) is positive, the pseudo-sound velocity can be used to distinguish between regions of quasi-neutrality and of charge separation in the expanding plasma. It is, hence, a keen indicator for the break down of the quasi-neutrality assumption. In section 6 we introduce a further formulation of the plasma expansion model, allowing the analytical investigation of the break down discussed in this section. We shall treat the special case of isothermality, y = 1, and charge separation (see fig. 1 and fig. 2). Before doing this, it is necessary to

Cli.

Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

349

6

a 5 ~t=15

3

2 ~t~9

t~O 0 -60

I

-50

-40

-30

I

-20

-10

I

0

.10

.20

Fig. 6. Square of the pseudo-sound velocity c as a function of x for three different time steps: (a) in comparison to the space-time evolution of the ion density n and the electric field E; (b) in the isothermal case, y = 1, with charge separation.

analyse the numerical results for this case in more detail. For this purpose the space-time behaviour of the solution near the collapse is shown on an enlarged scale in fig. 7. Figure 7a shows the gradual steepening of the ion velocity. If t approaches the collapse time = 17.94, the gradient of v becomes minus infinity. Figure 7b shows the electric field which also steepens and whose gradient goes to plus infinity for t —+ t,. The specific volume V = 1/n illustrated in fig. 7c yields some interesting details being of importance for the analytical investigation of the collapse. At any fixed time t the inverse density is V-shaped consisting of two branches with different slope which remains finite for t—* ti,. The minimum of V decreases linearly and becomes zero, i.e. n—* +so, at’ x = x~,= 4.112. At t = 17.935, e.g., the maximum value of the density equals 11.35. In the following we call the point (x~,t~)“critical point”. The electron density and the electrostatic potential which are not shown here, remain smooth and monotonically decreasing at the critical point. It is remarkable that the numerical procedure, consisting of thousands of iterations, is able to resolve these fine structures up to the breaking point. Furthermore, we mention that by using the (second)

350

Ch. Sack and H. Schamel, Plasma expansion into vacuum

2.6

4.00

4.04



A hydrodynamic approach

4.p8

x—’.

0.16

Fig. 7. Spatial dependence of the ion velocity v, (a), of the electric field E (b), and of the specific volume V~1/n (c) for four equidistant time steps just before wave breaking (see also fig. 1 and fig. 2); the minimum of the specific volume approaches zero linearly (broken line in (c)).

Poisson solver within the Lagrangian scheme, additional spatial derivatives which may become singular, are avoided (see section 4.2). In the (second) Poisson solver only the potential 4’ is differentiated, which is the smoothest quantity in the system. Therefore, the Lagrangian scheme together with Poisson solver II not containing a1n, is advantagous compared to the Eulerian scheme and Mason’s Poisson solver [38]. The structures we are treating here appear on a smaller scale than the Debye length. Thus the question may arise whether the basic model is valid globally. At the moment, however, we deal with the mathematical aspects of the dissipationless plasma expansion; its physical aspects will be discussed in section 7. The numerical results presented in detail in fig. 7 and the characteristic features deduced there from

av/ax~—so,

aE/ax~+so,

V~0 (n~so)

for

(x, t)~(x~, t~),

(5.5)

Cli. Sack and H. Scliamel, Plasma expansion into vacuum — A hydrodynamic approach

351

as well as the linear behaviour of V in the critical region are the starting point for the analytical investigation of the collapse. The solution shown in fig. 7c can be described in good approximation by the following empirical formulas: V(x, t) =

‘1(t)

t—

Xmin(t) =

~(t) t.

X(t)

[x



xmjn(t)]

+

. . .

(5.6)

with (5.7)

Equation (5.6) and (5.7) hold in the region x xmin(t), where xmjn represents the minimum of the V-profile. A similar ansatz with another .~K(t)can be made for the second branch x ~ xmjn(t). In (5.6) and (5.7) x and t represent transformed variables, corresponding to the original variables via the transformation x—*x~—x

and

t—*t~—t.

(5.8)

This transformation represents a reflection of the V-profile on the critical point (x~,t~)and a subsequent shift of the latter to the origin. The critical point is approached from the right, x, t—* 0~. Inserting the transformation (5.8) into Poisson’s equation (2.24) and into the equations of motion in the Lagrangian picture, eqs. (4.9)—(4.11), the transformation of the dependent variables reads v-*v,

V—*V,

4’—*4’,

E—*-E.

(5.9)

To lowest order in the smallness-parameter t [see (5.8)] the quantities ‘1(t), ~‘{(t)and ~‘(t) in (5.6) and (5.7) are non-vanishing constants. In the next section we shall justify the empirical formulas (5.6) and (5.7) and clarify the collapse mechanism.

6. Theory of bunching and wave breaking in ion dynamics The characteristic features of the numerical break-down [see (5.5)] described in the preceding section indicate the occurrence of wave breaking. The density becomes strongly peaked and finally diverges. As mentioned in the Introduction, this unbounded increase of the local ion density is called “bunching”. The relationship between this behaviour of the ion system and the dynamics of nonlinear Langmuir oscillations, which will be treated in appendix A, enables us to investigate analytically the collapse in the ion dynamics. In analogy to the analysis of nonlinear Langmuir oscillations (Davidson [62] and Kalman [73], see appendix A) we deduce in section 6.1 a single scalar wave equation. It describes the ion system including charge separation. Compared to the case of electrons, where one gets a simple oscillator equation for the electron velocity, this wave equation is a more complicated one and in general cannot be solved analytically. In section 6.2 a perturbative solution of this equation in the vicinity of the critical point (x~,t~)is presented for the isothermal case, y = 1. It reproduces all the properties of the numerical solution (see fig. 7) proving the existence of wave breaking.

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

352

6.1. Scalar wave equation In order to deduce the scalar wave equation for the ionic system, two important steps have to be performed. The first step consists of the transformation to the Lagrangian mass variable [see (A27) in appendix A]:

(6.1)

With reference to thermodynamics, in the second step the specific volume V = 1/n is introduced. It is useful to replace the density n by the specific volume, since in contrast to n, V is a limited function in the case of unlimited bunching (V—* 0 for n —~so). In differential formulation (6.1) reads: (6.2) where t is held fixed. With regard to the transformation made in section 5, eq. (5.8), is replaced by From (6.1) it follows that i~—*0’~’ for x—*0”. With (6.2) eqs. (4.9)—(4.11) transform to: —~

~.

x(~,t) = v(~,t)

(6.3) (6.4)

o(,~,t)

=

E(~,t)

=

-

V(71,

t)

~

4’(m

t).

(6.5)

For Poisson’s equation we obtain: [~

~-

~-

4’(~,t)]

= V(~, t) Tie(4’)

-1.

(6.6)

Equations (6.3)—(6.5) can be summarized to give

a~VasV=_a71[~ a~4’].

(6.7)

The addition of equations (6.6) and (6.7) yields: ne(4’)

=

(1-

(6.8)

~)/V.

Inserting for ne(4’) in eq. (6.8) the generalized equation of state, eq. (2.25), the potential expressed in terms of the specific volume V: -

~ y-1

[(1~ti1 t~ v i

i’

4’ can

be

69

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

where, at the moment, y ~ 1 is assumed. By means of (6.9), +

a~[~(1 ~)~2

a~

~)]

=

0,

4’ is

353

eliminated from (6.7):

1 ~ y ~2.

(6.10)

Equation (6.10) represents the scalar wave equation we are looking for. Its solution in connection with eqs. (6.3)—(6.5) yields the nt-dependence of x, v, E and 4’. Finally the solution in the Eulerian picture is obtained by inverting x(q, t)—* ~(x, t). In contrast to the electron dynamics where one obtains a simple oscillator equation [see, e.g. eq. (A30)], which can completely be solved analytically, the wave equation (6.10) for the ion dynamics is more intricate. It represents a nonlinear, partial differential equation in t and ~, which can be solved only by perturbative means. Moreover, the appearance of ~j-derivatives in (6.10) gives rise to a wave-like motion in the a-space so that in the stationary coordinate system always new particles, respectively fluid elements, start participating in the dynamics. This is, however, not the case for the nonlinear motion of electrons (Kalman [73]). There, a wave-like propagation along the 7)-coordinate is missing (a,1 = 0), and only those particles which are originally disturbed, participate in the bunching process. For y = 1 and y = 2, we get the following special versions of eq. (6.10): y=1:

v+a,1[~a,1(~/)]o

y=2:

V+

a,1[.~ a,1(

_V)]

(6.11)

=0.

(6.12)

Equation (6.11) is used in section 6.2 to investigate exemplarily the collapse in the ion dynamics by perturbation theory. In the quasi-neutral case, Tie(4’) from (6.8) must be replaced by ~e(4’) = 1 /V(asn). We obtain for 4’:

4’

y~i [(1)7]

By inserting

4’ from

(6.13)

(6.13) in (6.7) we get the quasi-neutral form of the scalar wave equation:

V+a~(1/V)~=0.

(6.14)

Even this much simpler scalar wave equation can generally not be solved analytically. Assuming especially a self-similar behaviour of the system, the solution of (6.14) can be written in a very simple form by means of the group theoretical formalism in section 3.1: V(~) =

(

/~2)1/(Y+1)

~as

~/t

where ~ is given by (6.1). In the isothermal case, y

V(fl=1/~~~, ~

(6.15)

,

=

1, the solution reads: (6.16)

354

Cli. Sack and H. Scliamel, Plasma expansion into vacuum — A hydrodynamic approach

It is a characteristic feature of the solution (6.16) that the second order t- and i~-derivativesin (6.14) vanish simultaneously. Because of V= 1/n we obtain for the ion density n( ~): n(fl=~~~, —1~~0,

(6.17)

i.e., in the c-space the ion density decreases linearly with the slope —1. Differentiating (6.17) via ~= ij/t = (1 /t) $~ n(x’, t) dx’ with respect to x we get

lan a 1 ——=—lnn=——. nax ax t

(6.18)

In accordance to (3.31) the x-integration of (6.18) yields the functional dependence of the ion density on the self-similar variable T = x/t. The next subsection is devoted to the perturbative solution of (6.11) in the vicinity of the critical point, i.e., for small values of ij and t. 6.2. Perturbative solution of the scalar wave equation for isothermal electrons We investigate the solubility of the scalar wave equation (6.10) by perturbation theory in the special case y = 1 [see (6.11)]. With some modifications the formalism used is applicable to other values of y, too. It is our aim to determine the functions ‘1(t), ~1{(t)and ~t) in (5.6) and (5.7), and to prove the behaviour of the dependent variables v, E and V, as it is shown in fig. 7. For the solution of (6.11) we make the following ansatz: V(~,t)

=

at[1



ht



(b

0



b1t)z

+

(c1

2

+

c2t)z

+

d 3 + 0(r4)] 1z

(6.19)

with z = lilt and 0< ij ~ lilt or lilt < z ~ 0, respectively; ~ = fit represents the position of the minimum of the V-profile in the ijt-space. The quantities (1, a, h, b 0, b1,. are constants. According to the transformation (5.8), eq. (6.19) represents the left-hand branch of the V-profile. A similar ansatz with other constants b1, c1,... holds true for the right-hand branch. We assume that t and z are quantities of the order r and determine the constants in (6.19) by comparing the coefficients in each power of r. It is shown that i) within each order of r the wave equation (6.11) is satisfied by the ansatz (6.19), ii) the solution represents exactly all properties of the numerical solution in the vicinity of the critical point (x~,t~), iii) to lowest order the solution is described by the equation —



.

a~v+ v a1v

=

.

0.

(6.20)

In the following, the limit t—*0, ij—~0,is made under the constraint ~q—* lit, i.e. we are at the minimum of V. A more general approach to the critical point, e.g., ‘q = KIlt, 0< 1, however, does not affect the essential results. At first it follows from the ansatz (6.19) that V—*0 for t—*0, guaranteeing automatically the divergence of the ion density at the critical point, ~ = t = 0 [see transformation (5.8)]. Since the ion K

S

Ch. Sack and H. Schamel, Plasma expansion into vacuum —

A

hydrodynamic approach

density ne remains finite at the critical point, it follows from (6.8), that for t—+ 0 and

355 ~j = flt—~ 0, V

must go to unity. In the ansatz (6.19) we perform the first and the second time differentiation: V’= a[1 + (b012 2h)t b0z b111t2 + 2(b1 c111)zt + c1z2 +(2c 2t 2c 2 + d 3 + 0(e4)] (6.21) 2 3d111)z 2Qzt 1z 12=a[2(b 2—/3 2+O(e3)], (6.22) 0uil—h)—a1t—a2z+/31t 2zt—/33z —











where a 1

=

2fl(2b1



(6.23a)

c1Q),

(6.23b) 2, 2c21’1 1~2 = 211(4c 2

(6.23c)



/33 = 2(3d1Q



3d111),

(6.23d)

c2).

(6.23e)

With V(~—+0~, t—*0~)—-*1we obtain from eq. (6.22) the coefficient h: h

=

Thus 1



b0ul



1/2a.

(6.24)

0(e), and it follows from (6.22) 2 + /3 2 + 0(e3)]. 1— V= a[a1t+ a2z— 131t 2zt+ /33z Furthermore, we get from the boundedness of ne the coefficient a V=

(6.25)

1 it holds [see (6.8)]: ne=(1~V)/V—+ne*;

t—+0~,

By means of (~.25)and (6.19) we get for a1

=

211(2b1



c1fl)

~=Qt—~~0~ ~,

.

(6.26)

t—+0:

(6.27)

fle*

ne* represents the finite value of the electron density at the critical point. In order to determine the other coefficients, we insert the ansatz (6.19) into the wave equation (6.11). By using (6.19) and (6.25) and expanding 1/V up to the order 0(1), we obtain explicitly the term: 2+ $ 2 + 0(r3)] [1 + ht + b 1 VV = a[a1t + a2z $1t 2zt + $3z 0z + 0(62)]. (6.28)

1



Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

356

For the ij-derivation of (6.28) one gets:

a,1 (1

[a2+g1t+g2z+0(r2)],

V)=t

(6.29)

where = /32 +

ha 2

g2

= 2(f33 +

+

(6.30a)

b0a1

b0a2).

(6.30b)

Next we determine from (6.25) the term 1/(1

L1

1 1—

=

1 [ a S(~,t)

S(ij, t)

=

a

with



V)



through Taylor expansion:

R(~jt) 21 S(~,t) + 0(r)]

(6.31)

0(e)

(6.32a)

2 + /3 2 = 0(62). R(-q, t) = —f31t 2zt + f33z The multiplication of (6.29) and (6.31) yields up to 0(1):

(6.32b)

1t + a2z

=

(6.33)

l_~,1(V)atS(fl,t)[a21t2z_a2S(~,t)+0(6)]

The n-derivative of (6.33) then yields the second term in the wave equation (6.11):

1

a,1[ 1 1

~

~)]

=

~

{—a~S + S {g

a,1(

2S +



a2 [(/32 + g1)t +

(2/33

+ g2)

z]} + 2a~R

0(r~)}.

(6.34)

For the sake of simplicity the arguments of S and R are dropped in (6.34) [compare (6.32a,Including b)]. 3 0(r~). Finally we insert (6.34) into the wave equation (6.11) and multiply with atS terms of order 0(62) we get: ‘-~

—cx~S+ S {g 2S



a2

[(1~2 +

g1)t+

(2f33

+ g2)

z}

+

2a~R + 0(r~) 0.

(6.35)

The demand that the wave equation (6.11) is satisfied to any order of r allows the determination of the coefficients. Comparing the coefficients with respect to r~and ~2 in (6.35), we obtain a2 = 0 and g2 = 0 from which by means of (6.23a—e), (6.27) and (6.30a,b) follows 2, b1

= ~

c1

=

~~eJ6Q

d 1

=

c2/3h2,

1~2 =

6c2u1.

(6.36)

For the evaluation of the term proportional to 0(r~)the order ~ in (6.19) must be taken into account,

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

357

which will not be pursued further. Therefore, we obtain by means of (6.36) and (6.24) for V(-q, t) in (6.19): V(~,t) =

at[1

+

~—



b~i~ + ~

(~2 —

112t2)

+

(~ lit)2 (ij + 211t) + —

~

(6.37)

To determine c, in (6.37) we use the fact that the derivative of the electron density 1 axne

= ~

1 anne

/1—1~\

)

~ a,1~\ ~

=

(6.38)

at the critical point is finite, too. Therefore, we obtain c 2

(6.39)

= ~bofle*/6hl•

The final solution V(’q, t) depends only on four constants a, b0, 11 and ne., which can be obtained from the numerical solution. To the lowest order, a represents the slope of the line V = Vmjn in fig. 7c; b0 denotes the finite gradient at the minimum of V. By means of (6.2)—(6.5) we now prove the singular behaviour of a~vand ~ at the critical point, = t = 0; from (6.37) we obtain to lowest order: 1

a~E=

1

1 ~

Vi 1

a~E=

~

1

..V

a,~x=

(6.40) (6.41)

—~+so.

-~ -~

Reversing the transformation (5.8) and (5.9), we get a1u—-+ —so, and ~,~E—++so in accordance with the behaviour shown in fig. 7a and fig. 7b. Equation (6.4), relating the Eulerian and the Lagrangian variables is used to transform from the ~jt-variablesto the xt-variables. The 7)-integration of (6.4) yields: x(~,t)

= Xmjn(t)



= Xmin(t) +

fd~’V(~’,t) at {(i

3

+

2

[~(n



+

f13t3)

(,~ —



fit)

112t2



~~0(2 —

Q2t2)

(~ lit)] + 0(r~)), —

0~

~ lit;

(6.42)

describes the time dependence of the minimum of V. By defining

Xmin(t)

~:=

[x



xmjn(t)]Iat,

(6.43)

Ch. Sack and H. Schamel. Plasma expansion into vacuum — A hydrodynamic approach

358

we obtain iteratively from (6.42) to each order the desired inverse transformation ~ = ij(x, t). This iteration procedure demands that ~ = 0(r) and x xmjn(t) = 0(62), respectively. The smallness of ~ means that the expansion is done in the vicinity of xmjn(t). The region of validity of the expansion, therefore, shrinks to zero for t—+0. Up to second order in e we obtain from (6,42) and (6.43): —



flt= ~[i +

(~

~ + ht)] +0(~~),

where h is given by (6.24). With (6.37): V(x, t)

=

at{i



~j

(6.44)

from (6.44) it follows for the x-dependence of the specific volume in

ht— b0 ~+ ~

(~+ 2ult)



b0

(~-~

+

ht)]

+

o(r3)}.

In a similar way we obtain by integration of the equations a~u= V and a~E= V with respect to replacing i~by eq. (6.44) the velocity and the electric field E depending on x and t:

v.

v(x, t)

= Vmjn(t) +

E(x, t) =

a~(1

Emin(t) + ~ +

+

t/2a)

+

0(~~)

(6.45) ij

and by

(6.46)

0(r~),

(6.47)

where Vmin(t) and Emin(t) are integration constants with respect to ~ or x. The latter determine the time behaviour of the velocity and of the electric field at the minimum of the V-profile, x = xmjn(t), corresponding to i~= lit in the n-space. By means of eq. (6.42) and of the relations ~ t) = v and ~(vj, t) = v(rj, t) = E(~, t), evaluated at the point i~= lilt, we calculate Xmjn(t) and vmjn(t): ~min(t)



alIt (1

t3min(t)



all

+



ht)

(6.48)

= Vmin(t)

lit (aboli



~+ an~)

(6.49)

= Emjn(t)

with Emin(t) = —E~+ 0(t), and Vm10(t) = v~.+ 0(t); v~,and E~ represent the finite values of the velocity and of the electric field at the critical point, respectively. In the two lowest orders the t-integration of eq. (6.49) resolved with respect to limin yields: 2), Vmin(0) = v~. 0(t From (6.50) we determine xmjn(t) by a further t-integration of (6.48): Vmin(t) =

v~+ (—E~+

(6.50)

all) t +

xmjn(t) = v~t+ ~(—E~ + 2ali) t2 + 0(t3),

xmjn(0)

=

0.

(6.51)

By comparison of (6.45) and (6.51) with the empirical formulas (5.6) and (5.7) we get ‘1(t)

=

a(1

~(t)

=

b



ht +...)

0[1 + (h



n~./3lib~) t +...]

(6.52) (6.53)

C/s. Sack ani”H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach =

v~.+ ~(—E~ + 2afl) t

359

(6.54)

with h from (6.24). Taking into account the third order in r for v(x, t) in (6.46), we get a term proportional to 2/a2t. Therefrom it follows that for t—+ 0 the second derivative of v with respect to x = [x xmin(t)] diverges, a2vlax2—+so, i.e., a2xlav2—+0. Hence, for the inverse function x(v) the critical point is a turning point. This result exactly corresponds to the characteristic behaviour of the velocity in the case of wave breaking of the hydrodynamic system introduced in section 3.2 [see (3.54)]. In a similar way we obtain for the electric field a2E/ax2 —+ so and a2x/aE2 —*0, respectively. From eqs. (6.46) and (6.47) and from the Eulerian equation of motion —

a~v+ v a~v=

E

(6.55)

we get a further proof for the existence of wave breaking. At first we verify that eq. (6.55) is satisfied up to 0(r). Inserting (6.50) into (6.46), we get x/t as the leading term, which is of order 0(1) because of x, t—0(e):

v(x, t) = x/t + 0(r).

(6.56)

To lowest order a~vand v a~vare 0(r~) such that eq. (6.55) reads

a

1v

+

v a~v= 0.

(6.57)

In the same way we show to lowest order by means of n(x, t)= 1/V(x, t)

=

1/at+ 0(1),

(6.58)

and

E(x, t) =

—E~ +

0(r),

(6.59)

that the equation of continuity and Poisson’s equation are satisfied. The leading term in (6.58) is 0(6_i) and reflects ion bunching. Equation (6.57) is the most simple version of an equation describing wave steepening and proves definitely that the ion dynamics is subject to wave breaking in the vicinity of the critical point. Just before wave breaking the electric field becomes unimportant, and the system behaves like an ordinary fluid. The dynamics of the system is thereby controlled by the convective term v a~v,which diverges. In a compressional wave, in which a region of high density propagates in the direction of low density, steepening and wave breaking are driven by a decreasing velocity profile. This is in accord with the numerical results shown in fig. 7 which were the starting point of our analytical investigation. Similar conclusions can be drawn from figs. 4 and 5. The case y = 1 and quasi-neutrality (see fig. 3) is an exception. Here, the velocity increases monotonically, and the conditions for wave breaking are not met. As mentioned in section 5, this, however, is due to the special choice of the initial conditions. If we had started in this case with a sufficiently decreasing initial velocity v(x, 0) or with another density profile, we also could have obtained wave breaking.

360

Cli. Sack and H. Scliamel, Plasma expansion into vacuum — A liydrodynamic approach

In the following consideration it is shown, that the special case y = 2 and quasi-neutrality (see fig. 4) is equivalent to other physical systems which are known to exhibit wave breaking. With tte(4’) = = 1 + 4’/2 [see eq. (2.21)], and the definition h(x, t) 2 n(x, t), the hydrodynamic ion equations (2.22) and (2.23) turn into

a~h+ a~(vh)= 0,

(6.60)

a~v+va~v+a~h=0.

(6.61)

In the ordinary theory of fluids the equations (6.60) and (6.61) describe the rotationless dynamics of long wave length, shallow water waves in a channel where h stands for the depth of the channel [55]. The solutions of these equations are the Riemann simple waves which, as known, show wave breaking (compare section 3.2). Therefore, it is not surprising that in our system wavebreaking exists, too, especially when we start with a more general set of ion equations. The equations (6.60) and (6.61) are, furthermore, identical to those describing the wave propagation perpendicular to a magnetic field in a small /3 plasma (/3 plasma (f3 = nT/(B2/8ir) with wave lengths exceeding the collisionless skindepth cIw,~[55]. In this case, the depth of the channel h has to be replaced by the magnetic field B and the sound speed by the fast magneto-acoustic speed. Summarizing, the collapse in the ion dynamics can be understood by the following simple picture. The ambipolar electric field, being finite at the front of the plasma and hump-like, accelerates the ions and generates a similarly shaped velocity profile. This profile continuously steepens due to the ion nonlinearity and, in a further stage, when the ion flow velocity becomes supersonic, the convective term in the ion momentum equation overrules the electric field. Accompanied by an unlimited density bunching this leads to a simple wave structure analogous to an ordinary hydrodynamic compressional wave. The main conclusion drawn from this section, therefore, is that the dissipationless expansion of a plasma into vacuum is subject to the phenomenon of wave breaking, a phenomenon which came into attention only recently. It is known from ordinary hydrodynamics that after the advent of wave breaking, the wave does not remain in a simple wave structure. In the configuration space the dependent quantities as e.g., the velocity and the density, are no longer uniquely defined [53]. The description of the plasma expansion within the framework of a dissipationless hydrodynamic model itself, breaks down. In order to investigate the ion dynamics further, there are essentially two concepts. In the first one the dissipationless regime is kept. This requires the transition to the microscopic Vlasov description. In the next section we shall discuss among others the kinetic structures arising within this description after wave breaking. It will be shown that in the case of strong nonlinearity our solution coincides with the kinetic one. In the second concept one keeps the hydrodynamic description and introduces a viscosity term in the ion momentum equation which prevents the unlimited steepening of the velocity profile and, therefore, the wave breaking. This procedure is well known in hydrodynamics [65, 66]. Our numerical investigations using the second concept are presented in sections 8, 9 and 10.

7. Comparison with previous hydrodynamic and kinetic models In the quasi-neutral case the breaking of ion acoustic waves with finite amplitude has been investigated by Akhiezer [32] and by Gurevich et al. [35,43]. Akhiezer’s model was based on the

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

361

kinetic description of the electrons using the Vlasov equation and on the hydrodynamic equations for cold ions. Gurevich et a!. used the ion-Vlasov equation at finite ion temperature and investigated their nonlinear self-similar behaviour; thereby the electrons were assumed to be isothermal. Ignoring for a moment the fact that quasi-neutrality is violated at the steepened front, we can get analogous results by making use of eq. (6.14). By solving eq. (6.11) explicitly we have shown that the structure of a simple wave does not rely on the assumption of quasi-neutrality which loses its validity prior to wave breaking. Simple waves can be found in the case of charge separation, too. In some publications [34, 35, 43], it was conjectured that the inclusion of dispersion, which in the description of ion waves comes into play through the term a~4’~ 0, can prevent wave breaking. Our results do not support this conjecture in the extremely nonlinear case of plasma expansion into vacuum. In addition the collapse cannot be prevented either by finite ion temperature in the realm of a collisionless model [56]. From the numerical point of view the resolution of wave breaking in connection with ion bunching is only possible within a numerically dissipationless discretization procedure, e.g., the explicit Lagrangian scheme [65,66]. Explicit algorithms fixed in space such as Lax schemes are less suitable. The associated numerical dispersion and dissipation [66] cause an unacceptable strong smoothing of the nonlinear structures. Such a procedure was used by Widner et al. [37], who were the first to investigate numerically the isothermal hydrodynamic plasma expansion into vacuum, taking into account charge separation (see fig. 8a). The limitation of the ion front velocity to about 3c~found by them together with the ion peak, is an indication for the onset of wave breaking. However, as mentioned by Crow et a!. [42], the influence of inaccurate boundary conditions for Poisson’s equation (4’ = const. at the boundaries) cannot be excluded. An improved description of the plasma expansion with respect to the boundary conditions was given by Crow et al. [42].The hydrodynamic ion equations were solved by the explicit Lagrangian procedure similar to our method. The main difference to the expansion model discussed in this review is the choice of the initial ion density profile which was a step function. Thus 4” is discontinuous at the step which survives for all times. Using the continuity of the potential and of the electric field at the front, the solution in the electron cloud region is matched to the numerically obtained solution in the ion region. In the early stage of the expansion one observes a well-pronounced hump in the ion density which disappears for large times. The density profiles in this stage are similar to those shown in fig. 1. Crow et al., however, do not make any statements about the appearance and disappearance of the hump. The space-time history of their solution for the ion density is illustrated in fig. 8b. It is possible that the disappearance of the ion peak in their solution is due to the fitting procedure. From hydrodynamics it is known that the application of jump and continuity conditions, respectively, in the region of discontinuities introduces dissipation [53]. This smoothes small scale structures and reduces amplitudes. Shock fitting, in general, has some undeniable shortcomings [65]: Firstly, shock fitting can only be used if the solution is known exactly in front of the discontinuity as in the model of Crow et a!. Usually this is not the case if in an initially diffuse profile discontinuities appear [65]. In gases or fluids, shocks can develop spontaneously, and there is no way to find out when and where shock fitting should be applied. Secondly, the finite grid size does not allow to localize the shock front precisely. This introduces inaccuracies of the dependent variables, falsifying the numerical solution. Crow et a!. do not specify how they have tried to circumvent this problem. In a more recent publication, Gurevich and Meshcherkin [49] working out the model of Crow et a!. once more, presented identical solutions. Their ion hump shown in fig. 8c, was not commented. A similar model was used by True et al. [50]. The ion hydrodynamics was solved in the Lagrangian description by means of the conservative leap-frog method [66]. The stability analysis of the latter method leads to two uncoupled difference meshes which may drift apart in the course of time [66]. Such

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

362

n,n

a

I

—80

::

o.:.

~

—40

0

.0,1-I

40

80

120

x

180

0

50

100

x

150

~n

c

0.5

—1

0

0 1.0

—2

—50

1

2

3

4

r=

d

00

—~—

cot

0

50

100

150

Fig. 8. Space-time behaviour of ion(electron) density profiles taken from: (a) ref. [37],(b) ref. [421,(c) ref. of density profiles a, b and c with those shown in Figs. 1 and 2.

200

250

300 x 350

[491 and (d) ref. [50].Note the similarity

a decoupling can be avoided either by defining the dependent variables on one net only or by introducing diffusion that couples both meshes. The paper of True et al. does not explain which possibility was chosen. The striking smoothness of their density profile shown in fig. 8d, one may suspect, indicates strong numerical dissipation. With respect to the ion hump found by Crow et a!. [42], True et al. simply state: “We have not observed such a bump with the given initial conditions”. It is remarkable that four articles solving the same set of equations give essentially three different answers. As discussed we attribute these discrepancies to the uncontrolled action of dissipation. To study nonlinear structures within dissipationless ion hydrodynamics, a Lagrangian dissipationless scheme is preferable, like the one presented and reviewed in this paper. A scheme like ours offers the correct numerical answer to the plasma expansion problem, at least for diffusive step-like initial profiles. Finally, we discuss in this section the phenomenon of wave breaking within the framework of the collisionless microscopic model. Using particle simulations, the breaking of finite amplitude ion waves

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

363

has been investigated by Forslund et al. [74], and by Buchelnikova and Matochkin [75].In both works a homogeneous plasma is modulated initially by a sinusoidal ve!ocity perturbation with contro!able strength v0. Both simulations show- the occurrence of wave breaking in connection with ion bunching. As a result, kinetic structures arise, the properties of which are depending on v0, T1ITe, and the electron polytropic exponent y [74,75]. For small values of v0 (v0 ~ c,), double layer structures of the slow ion acoustic type [76, 77] are seen. In this case, the region between them is characterized by a laminar increase of the ion temperature [78]. For larger values of v0 (v0 ~2—5c~)a symmetric x-point structure in the ion phase space emerges. For still larger values,of v0(v0 5—10 c5) a mixed type of wave breaking structure is obtained, in which ion reflection and free streaming occurs simultaneously (see ref. [74] fig. 7 and ref. [75] fig. 18). For velocities far beyond ion acoustic speed (v0 ~ 20 c5) an inverted S-type structure arises expressing multiple ion free streaming. Figure 9, taken from Forsiund et a!. [74], illustrates this specific pattern. In this case charge separation limits the electrostatic potential 4’, so that the conditions for ion trapping and reflection cannot be reached. In a simple analytical model, Forslund

:~* I

200

~

Fig. 9. Ion wave behaviour taken from ref. [74]for v~= 30 c, and TI T, = i0~ (a) ion phase space v versus x; (b) ion density; (c) electrostatic potential; (d) electron density; (e) ion phase space v versus x after wave breaking. The wave has broken by multiple ion free streaming.

Ch. Sack and H. Schamel. Plasma expansion into vacuum — A hydrodynamic approach

364

et a!. show that the electrons cannot effectively shield the ions over distances smaller than several Debye lengths. The electron density lags behing the ion density. As long as the ion velocity is unique, the curves of n, v and 4’ obtained in these particle simulations, coincide with our hydrodynamic results. By comparison with the kinetic results, we can decide which kind of phase space structure would be approached after wave breaking if a kinetic treatment had been used. In the quasi-neutral case the potential grows with n, and ion reflection will be the saturation mechanism. Therefore, we expect in a kinetic description the transition to two double layers or to a x-point structure which may be modified by free streaming. If charge separation is taken into account, 4’ and ~e(4’) are more moderate and do not follow the increase of the ion density. Consequently, ion free streaming with triple velocity will come up after breaking, resembling the case of ultrasonic v0. The expansion of a plasma into vacuum can be interpreted as the extreme case of a density modulated plasma in which the density at the minimum reaches zero. In this sense, the plasma—vacuum system belongs to the strongest nonlinear systems one can imagine. Ion bunching and ion wave collapse are natural events in such a system. Appropriate kinetic calculations of the plasma—vacuum system are not yet available. Particle simulations are not suited if steepening occurs in the !ow density region (y 1) because of the enlarged numerical noise at low densities. Codes, solving Vlasov equation by direct discretization do not contain such kind of noise and seem to be more advantageous to investigate nonlinear phenomena associated with plasma expansion into vacuum. One has, however, to be sure that numerical dissipative effects are minimized by the choice of the difference method. One should not refer to LAX schemes [66, page 191 ff.]. In any case, the entropy production f dx dv f In f has to be calculated to check the rate of dissipation being aware that dissipation cannot be eradicated completely. In the hydrodynamic description, one way to get rid of the collapse phenomenon in a controlable manner, is to introduce a Navier—Stockes viscosity. Its physical foundation and its impact on the numerical procedure will be discussed in the next section.

8. Navier—Stokes viscosity and the implicit Lagrangian scheme The particle simulations consulted in section 7 indicate a thermalization of the ion distribution after steepening, pointing towards micro-instabilities. These anomalous collision effects are modelled by adding a Navier—Stokes term to the ion momentum equation (2.23), (4.11) respectively: a,v

+

~ a~v dv/dt = —a~~ + a~v,

= const.

(8.1)

being aware of the shortcomings of this approximation. A further justification results from the fact that the energy of the directed plasma flow is transferred to smaller scale-length bringing viscous effects into play [55].This assumes that T. is not completely negligible [79]. The numerical results in sections 9 and 10 will show that the behaviour of the expanding plasma is rather insensitive to the strength of v. Introducing the Navier—Stokes viscosity, the ion momentum equation becomes parabolic in space allowing the implicit formulation of the discretized set of equations (see, e.g. [66]). Accordingly, the

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

discretization of

p

2

365

a~vin (8.1) reads

2 n+1/2

2

2 n—i/2

(1 a) (~v/ax ). ], (8.2) involving both, the new and the old time step. In (8.2) j and n represent discrete space and time !abels; n+1/2 a denotes the implicitness parameter, 0 a ~ 1. The explicit scheme is given by a = 0, in which is determined a!gebraica!ly by the known quantities at the previous time step. For a 0 the actual time is introduced, and the scheme is called implicit. The special, choice a = ~ yields the half-implicit or Crank—Nicholson scheme [66]. Of special importance is the full-implicit scheme, a = 1, which is used to get the results presented in sections 9 and 10. With a 0, eq. (8.1) is transformed to a set of linear equations for v’~~2 characterized by a tn-diagonal coefficient matrix. The first and the last row of the latter have to be completed by the boundary conditions. For this purpose we use the discretized version of eqs. (4.33) and (4.44), respectively from section 4.3. Further details of the implicit difference scheme are discussed in appendix C. An essential point in the implementation of evolution equations by finite differences is the stability property of the algorithm. Such an algorithm is stable if small computational errors do not accumulate in the course of time [66]. The stability property can easily be established for linear systems of partial differential equations. Basic difficulties arise in non-linear systems. In hydrodynamics or magnetohydrodynamics, e.g., the pressure gradient term —(1/n) a~pis replaced by —(c2In) a~nwhere c(x, t) dp/dn~,denotes the characteristic sound velocity. The closure of the system is achieved by assuming a polytropic equation of state p n7 corresponding to isentropy, s ln p/n7 = const. Keeping c2/n constant with respect to x and t, the pressure term becomes linear allowing a stability analysis (see [65,66]). Our system is similar. The electrostatic potential 4’, replacing p, is known by solving an additional equation (Poisson’s equation) and can be considered as a functional of the ion density n. Based on this concept we have introduced in section 3.2 the pseudo-sound velocity c, eq. (3.46), which in the quasi-neutral case becomes c = V5~~ In the homogeneous plasma (n 1) c corresponds to the usual ion acoustic velocity, and the linearization method is applicable here, too. In the charge separation case, however, c2 can become negative and c imaginary (see section 5). There exist obviously regions where the concept of pseudo-sound velocity and the related stability analysis lose their validity. In which sense the stability conditions describe the behaviour of our numerical model will be discussed next. The ion momentum equation (8.1) is first formulated in the Lagrangian coordinate ~: v[a(a v/ax ).

+





dv/dt =



a~V+ A a~v+ B a~v,

C

(8.3)

where the ~- and t-dependent coefficients are defined by c2V

2

0 V V V C:=—~-, A:=V-~af-~ and B:=p—~. In (8.4) c2 represents the pseudo-sound velocity in terms of ~ c2

(8.4)

= EV2/V 0

a~V

(8.5)

Cli. Sack and H. Scliamel. Plasma expansion into vacuum — A hydrodynarnic approach

366

Linearization of (8.3) is achieved by assuming locally constant coefficients (see, ref. [65, page 297]). The discretized version of eqs. (4.9), (4.10) and (8.3) gives rise to a closed set of linearized equations, the stability of which can now be studied. Details of the stability analysis are treated in appendix D. Here we discuss the results only. According to the sign of c2, two cases have to be distinguished. First we consider the general stability condition, inequality (D18) in appendix D, for the case c2 > 0, and three different values of the implicitness parameter a (see (Di9)). For a = one gets p2sin2kh~4

resp.

p2~4,

(8.6)

if one chooses the maximum value of the sinus (kh = ir/2) where p2 C2(V~ 2and h as ~ Since 5 .~t/Vh) in this half-implicit case (8.6) is independent of ~‘, the viscosity does not contribute to the stability of the scheme. This property can be verified numerically. In fig. 10 the ion front is plotted as a function of x at different times for v = 1. Ion bunching is weak, but in the front region oscillations from grid point to grid point develop increasing in amplitude. Shortly after this time, the solution breaks down. Obviously this break-down is of numerical origin in contrast to the collapse described in sections 5 and 6. Inserting typical values of the front region, the stability condition (8.6) is violated. We note that inequality (8.6) coincides with the stability condition for the explicit, inviscid Lagrangian scheme. The case a = 0 (a = 1) can be treated commonly (see (D19)): ~)

[1 (a

fl2

+ (fl2]

~ 1,

(8.7)

where the lower (upper) sign holds for a = 0 (a = 1). In (8.7) the stabilizing influence of v for the full-implicit scheme (a = 1) is visible. In the explicit scheme (a = 0), on the other hand, viscosity deteriorates the stability. Within the explicit scheme, oscillations like those for a = 0.5 (see fig. 10) arise. The solution breaks down the earlier, the larger i.’ is. Generally speaking, inequality (8.7) is in agreement with the behaviour of our code in case where the density decreases with x (c2 >0). 0.12

Fig. 10. Spatial dependence of the ion density for four different times obtained from the semi-implicit code, a separation. The fluctuations in the front region at t 21.9 are due to a numerical instability.

=

0.5, v = 1, y

=

I and charge

Ch. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

367

In regions of increasing density, i.e. c2 <0, one obtains an expression which formally implies that the numerical scheme is unconditionally unstable (see inequality (D20) in appendix D). This is accomplished easiest if a = 1. Accordingly, the full-implicit scheme should become unstable immediately after ion-hump formation in the region of increasing density. This is, however, not in accordance with our numerical experience. Our essential conclusion is that the stability analysis, based on the concept of pseudo-sound velocity, yields reasonable results in those regions where c is a real quantity and, thus, can be interpreted as a velocity. It must, however, fail if c is imaginary and consequently loses its physical relevance. It is indeed questionable whether stability conditions derived in such a limited sense, describe the actual behaviour of a numerical model for an extremely nonlinear and inhomogeneous problem. Strictly speaking, the set of discretized equations is a multi-dimensional, highly nonlinear discrete mapping. Its properties should be derived within the framework of the methods for investigating nonlinear dynamical systems (compare, e.g. [80]). Especially, it is contradictory to keep the quantities A, B, C, in (8.4) constant if, at the same time, the a-dependence of af V and afv is taken into account. Such an inconsistency holds for any linearization procedure of nonlinear equations. Conventiona! stability methods including hydrodynamic ones can at best yield a rough estimate of the dynamical behaviour; they become the more unreliable the more nonlinear and inhomogeneous the problem is. By means of the full-implicit Lagrangian scheme, a = 1, one is able to prevent wave breaking and to follow the expanding plasma over long times without any numerical instability. In the next section the numerical results for the early stage of viscid plasma expansion are presented paying special attention to the dependence of the solution on the viscosity v, on the charge conditions (charge separation, -neutrality) and on the initial conditions.

9. The early stage of viscid plasma expansion Figure ha shows the spatial dependence of the ion density in the isothermal case (y = 1) with charge separation and with u = 5 for three different time steps. In fig. lib the density n, the ve!ocity v, and the electric field E in the ion front region are plotted on an enlarged scale for t = 60. Different to the

__

a

-200

-100

0

100

x

200

Fig. ii. Ion density as a function of x for (a) three different time steps, for v = 5, y = 1 and charge separation; (b) shows the spatial behaviour of n, E, v at the leading front for t = 60. The dots representing n are the actual grid points.

368

Cli. Sack and H. Scliamel. Plasma expansion into vacuum — A hydrodynamic approach

corresponding inviscid case (see figs. 1, 2 and 7) where the solution broke down at t~,= 17.94, wave breaking does not occur here. The ion hump, which expresses the density bunching, increases in the course of time, but the rate of increase is not explosive. Thus, the implementation of a controllable viscosity takes care of stabilizing the ion bunching process. Within the full-implicit discretization scheme, a = 1, no numerical instability develops, although the ion hump is very spiky and steep as shown in fig. lhb. Each point in fig. hib corresponds to a grid point, i.e. the ion hump is resolved by a large number of grid points. For higher values of v, e.g. ii = 10, the ion bunching is further delayed (compare ref. [57]). The qualitative structure of the ion front appears to be independent of the numerical value of v. The influence of the chosen charge model on the viscid solution is depicted in fig. 12 for the case of y = 1.2, a = 1 and x = 5. Figure 12a shows the ion density at t = 60 resulting from charge separation (solid line) and from charge neutrality (broken line). The density curves differ significantly in the front region whereas in the rear they are indistinguishable. Figure 12a reflects the delaying influence of charge separation on the ion front velocity, confirming suggestions in this direction, which have been made by other authors too (see e.g. Denavit [46]). Behind the front the plasma behaves quasi-neutral. In this region the ion motion can be described by the self-similar theory (see section 3.1). 1.0’

0.75

\

—200

a

-100

0

100

~__._.

b 0.3

0.3 25

/

V

~

/

2.35

~1’ ... ‘~

2.3 84

~\~

240.2

v 0.1

-2

‘~~j \

0.224 n

—2.2

~

n

2.6

V

20

-4

0.1

~.

~

88

90

.

107

109

111

113

115

Fig. 12. (a) Ion density as a function of x for y = 1.2, v = 5 at t = 60, comparing the effect of charge separation (solid line) with that of charge neutrality (dashed line). (b) and (c) show the. spatial dependence of n, cb, v at the front of a larger scale.

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

369

A further comparison of the expansion pattern with and without charge separation is presented in figs. 12b and 12c, where n, v and 4’ in the front region are plotted as functions of x for t = 60. In the case of charge separation, fig. h2b, the ion density is needle-like, whereas the velocity v and the potential 4’ remain relatively smooth. Consistent with the charge separation model the width 4 of the Such a consistency ion front is comparable to the local Debye length given by AD(x, t) = t)~~2~2. is missing in the quasi-neutral case, fig. 12c, where the front width does not satisfy the necessary condition 4 ~‘ AD(x, t). In contrast to fig. 12b no pronounced ion hump appears. Instead a sharp density step is formed. Disregarding that the quasi-neutrality assumption breaks down at the density step, the behaviour of the dependent quantities in fig. h2c resembles a usual hydrodynamic shock. This is, however, not surprising, because in the quasi-neutral case the set of equations coincides with that of an ideal neutral gas (see section 3). In order to investigate the dependence of the solution on the initial conditions, two calculations have been performed, in which the initial velocity of the ions was unequal zero. The corresponding density profiles in comparison to v (x, 0) = 0 (broken line) are given in fig. 13 for y = 1.2, i.’ = 5, and charge separation at t = 60. The solid curve results from the initial velocity fle(X,

v(x, 0)

=

v 0[1



n(x, 0)] ,

(9.1)

where v0 = 1; n(x, 0) is the initial density profile from eq. (2.28). The dotted curve shows a density profile following from v(x, 0)

=

\/—2 4’(x, 0);

(9.2)

4’(x,0) is the self-consistent potential calculated at t=0. Both expressions (9.1) and (9.2) try to describe a plasma state in which the diffuse initial density profile, due to some complicated and in detail unknown generation mechanisms, is correlated with an initial momentum of the particles. Figure 13 clearly reflects the dependence of the front behaviour on the initial conditions. The number and the present position of the fast ions essentially depend on the details of the initial conditions. An important quantity for the description of the dynamics of expanding plasmas being of interest, especially in experiments, is the velocity of the fastest ions at a given stage. Figure 14 presents the time

-200 Fig. i3. Ion density as a function of x at line); u(x, 0) = V—2 cb(x, 0) (dotted line),

t(y =

60 for three choices of the initial ion velocity = 1.2, = 5, charge separation).

p

(t

x =

0): v(x, 0) = 0 (dashed line); v0(1 — n0(x, 0)) (solid

370

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

8

7

6

5

4

~

______

9

0

I

~

0

10

20

~)

40

50

60

70

90

Fig.14. Maximumionvelocityasafunctionoftfor1:y1,QN,v=0;2:y’~1.2,CS,i’5,v(x,0)V—2~(x,0);3:yi,CS,v0;4:yi. CS, p1;5: y=i,CS, v=5;6: y=1, CS, v10;7: y=1.2, CS, r5, v(x,0)=v 0(1—n(x,0)), u0—1; 8: yi.2, QN, v=5;9: y1.2,C5, v = 5. QN and CS denote quasi-neutrality and charge separation, respectively; except for curves 1 and 3 a = i is used in the calculations (_* full-implicit scheme).

evolution of this maximum ion velocity for a series of different cases. (Notice that vmax should not be mixed up with the final velocity of the fast ions, which will be investigated in subsection 10.2.) Firstly, it is striking that the expansion elapses mainly in two phases. The first phase is characterized by a continuous increase of Vmax~This increase results from the initial acceleration of the plasma, which is particularly strong, if the plasma was at rest at t = 0, i.e. v(x, 0) 0. If, due to the initial conditions, the electric field is hump-like, this behaviour transfers to the velocity profile as well in that first phase (see curves 3—6, 8 and 9). At the time where the solution of the dissipationless model breaks down due to the unlimited ion bunching, and the following wave breaking, the maximum ion velocity reaches in the dissipative case, v ~0, a constant value of several c5, depending on i.’, on y, and on the charge conditions. In this second phase the ion bunching is stabilized by the viscosity. For the sake of comparison curve 3 shows the maximum ion velocity in the inviscid case for ‘y = 1 with charge separation. The cross at the end of the curve marks the critical time t~,= 17.94, at which the solution breaks down (see sections 5 and 6). If the plasma has an initial velocity at t = 0, the acceleration is diminished due to the reduced electric field (see curves 2 and 7). In the case where the initial velocity is not adapted self-consistently, eq. (9.1), a time evolution of the velocity similar to the cases with v(x, 0) 0 can be seen (see curve 7). Different to this, the use of a self-consistent initial velocity profile, eq. (9.2), leads to a weak, but Vmax(t)

C/s. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

371

continuous increase of the maximum ion velocity (see curve 2). Therefore, no density bunching takes place in the considered time interval (see also dotted line in fig. 11). In fig. 14 a special case is the quasi-neutral one with y = 1, v = 0 and v(x, 0) 0 which shows a strong continuous increase of the maximum ion velocity (see curve 1). This behaviour can be understood by the following argumentation. Using quasi-neutrality, ne = n, and an electron equation of state with a polytropic exponent y> 1, one obtains from eq. (2.25) the relation ~E=a~(n71).

(9.3)

For x—* +so the right-hand side of eq. (9.3) vanishes because n and a~nboth go to zero. As long as y ~ 1, the electric field also vanishes. A finite value of E in the vacuum region, which is necessary for a strong continuous acceleration of the plasma, can only be achieved for y = 1. In this special case it holds E = —a~ln n, and a constant value of E is necessarily correlated with an exponentially decreasing density (see also fig. 3 and section 5). This is exactly the case here. In all the other cases the electric field vanishes for large x within our simplified model, and consequently, the conditions for ion bunching and wave breaking are met. In this section the acceleration and the bunching phase and its dependence on various parameters have been investigated within the dissipative model. An essential result is the stabilization of the ion bunching by the dissipative term being introduced in form of the Navier—Stokes viscosity. Furthermore, it has been found that the stabilization of the bunching is correlated with a reduced maximum ion velocity. One could believe that this is the final state of the expanding plasma. In the next section it is shown that the two phases discussed here are followed by a third phase in which the plasma is accelerated again. A simple explanation is given for the transition of this third phase. For long times the ion front velocity approaches the value given by the self-similar theory (see eq. (3.29)). Despite the fact that oscillations of the dependent quantities emerge in the rear of the ion front, the main part of the plasma behaves quasi-neutral and can be time-asymptotically described by the self-similar solution. The existence of these so-called intermediate asymptotics is proven by a perturbative analysis of the single scalar wave equation for the isothermal quasi-neutral case (see eq. (6.14)) using the self-similar solution as a starting point.

10. The long-time behaviour 10.1. The debunching process

To get an impression of the late phase of plasma expansion we made several long-term calculations which go beyond the bunching phase. Figure 15 exemplarily shows the evolution of the ion density in the time interval 0~t~140 for y = 1, v = 5 and 1= 1 (see eq. (2.28)) at nine time steps. One realizes the stabilizing influence of viscosity. The existence of the finite peak and the associated fast ions is the relic of the nonlinear process taking place in the early stage. At t = 95 (see dashed line in fig. 15) the peak reaches maximum height. Shortly after this the peak decays rapidly, the decay time being shorter than the time for the formation of the peak. The mechanism responsible for this decay can be understood by mass conservation and by the characteristic change of the velocity profile in the front region. This is illustrated in fig. 16 where the x-dependence of the ion density n, of the velocity v, and

372

Cli. Sack and H. Schamel. Plasma expansion into vacuum — A hycirodynarnic approach 14Q~ 120 100 40 20

\

1.00

\

\ 118

\

/

\

0.6

\ 0.4 —

0.2 0 -200

I

I

0

2(X)

Fig. 15. Space-time evolution of the ion density for v whereas for t >95 debunehing takes place.

=

I

1

400

600

—.x

1 in the time interval 0 ~ t~140; at r

S. y

=

95 (dashed curve) ion’ bunching is maximum,

as

b 0.4 n

-A

E v 0.3 .3

,~

a

~

0.1

~

268~

—r— 0 256 260

.3

0.2

.i

as

-0



,~‘

254268

272

V

-0.1 1.90

~

C

0,4

,--

V

A

0.4

E 03 3

03

E

Q~ -2

0.2

\.~

o.i

1

at

0

‘0

274



284

A

~

C

.

292

-0.1 1.95

0.2 2

at

288

v

i

296 -0.1

1.100

Fig. 16. Spatial dependence of the ion density n (solid line), of the electric field E (dashed line), and of the ion velocity v (dotted line) in the front region for (a) t = 90, (b) t 95, and (c) t = 100; at t = 95 the velocity maximum overtakes the density hump initiating thereby the debunching process.

of the electric field E is drawn on an enlarged scale for the time steps t = 90 (fig. 16a), t = 95 (fig. 16b), and t = 100 (fig. 16c). For t = 90, the velocity maximum Umax lies behind the density peak, causing a narrowing of the mesh point distance 4 in the ion front region. Due to mass conservation n4 const. (sections 4.1 and 5), the ion density is locally enhanced but the collapse is prevented by viscous dissipation. At t = 95 the velocity maximum overtakes the ion peak, being in front of it at t = 100 (see figs. 16b and c). This shift of Umax can be understood as the desire of the plasma to approach the

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A liydrodynamic approach

373

100

10

0

~

0

I 10



v

o

~

=

00~i~ 8

I

I.

I

0.2 1

0

Vf

I

p

oy5

i

6

0

V

°V



.~

ovrl

I

~.

V

V

V

0.2 ~

V. I -2 a 10 1

It

It

I

* 111114$

5

• °

1

10

20

50

100

t

.~

,~—

200 300

2

=

1; •: v

=

I

0

I

b

0

I t.

Ito

.....

1

.

5

Fig. 17. Time dependence of (a) the maximum electric field Ema.~and (b) of the ion front velocity, v for three different values of v: x: v = 5, 0: r’

••~• .~‘‘~‘ 6~

10

20

50

100 ~200 300

1, on a double-, resp. single-logarithmic scale 0.2; the three plasma expansion regimes are separated by t~and

t0.

self-similar state, being characterized by a linearly increasing velocity profile (see eq. (3.32)). Therefore, it is not dissipation which causes the sudden decrease of the density peak but a process called “debunching”. In the latter, the bunching is reversed, the mesh point distances being rectified. Consequently 4 increases and n decreases. After the passage of Umax through the ion front, the curvature of &‘, being always negative in this region, diminishes and the electric field can act more effectively on the dynamics. This instant marks the onset of a new, as it turns out final, acceleration phase. In fig. 17a the time behaviour of the maximum electric field Emax is plotted on a double-logarithmic scale for three different values of the viscosity v. Figure 17b shows the ion front velocity v~which is the velocity at the position of Emax (i.e. at the position x~where ~e equals n) on a single-logarithmic scale. Three phases can be distinguished. The first one is given by the initial acceleration of the plasma resulting in steep gradients of v and E, associated with ion bunching. At t = t~,,where in the dissipationless case the collapse would occur, the second phase starts. It is characterized by an almost constant front velocity Vf of about 3 c~(see fig. 14 for more details) and by a nearly constant electric field Emax. During this phase the ion peak is present. Debunching sets in at t = t~where the maximum ion velocity overtakes the density peak. Therefore, debunching terminates the second phase called for short “plateau regime”. The last phase is determined by the decrease of the maximum electric field and by the continuous increase of the ion front velocity. Small differences in the behaviour of Emax and are visible for different v. With increasing v, debunching is slowed down, and the electric field is slightly enhanced. Nevertheless the velocity Ut 15 reduced a little bit owing to the fact that the acceleration is determined by E + v a~vsuch that the increase of E is compensated by the negative viscosity term. The larger electric field gives rise to a stronger increase of Uf. Later on the influence of t a~Ubecomes negligible because the curvature of v disappears. This is true independent of the polytropic exponent y. fl,

Uf

10.2.

The asymptotic velocity offast ions

Figure 18, taken from ref. [81],shows on a semi-logarithmic plot the time-dependence of the fast ion velocity Vf for various values of the electron polytropic index y (ii = 5). The scenario of ion accelaration

374

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach



v~*(12)

10-

.: 2

0-

,,,,,,,,,

1

5

10

,,,l,,,,,,l,,

20

50

v~(1,5) v,~’(1.7)

~

-,

~

VV

(18)

Fig. 18. Fast ion velocity c 4 as a function of time ton a single-logarithmic scale for several values of the polytropic exponent y; the approach of v4 to the corresponding self-similar values u’ is evident.

presented in section 10.1 for y = 1 is negotiable for y > 1, too. The most striking feature of this figure is the approach of v~to a constant velocity for ~ x• This constant velocity turns out to be identical with the front velocity deduced from self-similar theory. The latter has been derived in section 3, eqs. (3.24)—(3.30). From the inequality (3.39) it follows that the self-similar front speed is given by 2i/~/(y 1) (10.1) —

.

U~(y) =

It is exactly this speed which is asymptotically approached by the front of the expanding plasma. With other words, for long times charge separation effects at the front retreat and the deviations to the self-similar state becomes negligible. The approach to the self-similar state will be outlined further in subsection 10.3. For ‘y = 1 there is no limitation of the front velocity, a case which seems to be in accordance with experiments (see section 11). For y > 1 the driving force of electrons is reduced and the front acquires a finite velocity. As mentioned in ref. [81] the knowledge of this final ion front velocity and o-f the front velocity during the plateau phase yields an information about the thermodynamic behaviour of electrons during the expansion process expressed by y and about the initial electron temperature Teo in the unperturbed plasma. 10.3. Intermediate asymptotics This approach to the self-similar state sets in at an early stage of the evolution in the rear of the plasma as mentioned already in section 5 (figs. 2, 5, 6). Figure 19 presents the space-time behaviour of the potential 4’, corresponding to fig. 15 (y = 1, v = 5). All 4’ curves have an exact fix point, 4’ = —1, at x = —40, the position of the front at t = 0. The linear part in 4’ broadens in the course of time. It represents the self-similar state given by (10.2) and extends almost up to xf, the position of the ion front. At the front the curvature of

4’

is

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

1~

~200

2~0

375

1~

490

6~

2

Fig. 19. Electrostatic potential as a function of x for five different time steps. The positions of the ion front x

4 are indicated by arrows. For late

times, the bulk plasma behaves quasi-neutral (straight lines of the potential curves).

non-negligible and changes from a positive sign at the beginning to a negative sign at later times. Negative curvature corresponds to a surplus of ions indicating the presence of the ion peak. After debunching which takes place at about t 95, an oscillatory behaviour of 4’ appears behind the front. In contrast to the small scale length at the front in the plateau regime (t ~ 95), the scale lengths of the oscillations behind the front have effectively widened. This is also seen in fig. 20 where both the self-similar quantities, labelled by the index ss and given by (3.24)—(3.27), and the actual numerical quantities are plotted as functions of x for t = 175. Note the characteristic deviations between them. The electric field has a small scale structure at the front where n is appreciably different from ~~e’ and exhibits oscillations of longer wave-lengths in the rear where n

600

700

~

800

t

Fig. 20. Ion density n, electron density n~(dotted line), electric field E and ion velocity v for v = 5 at = 175; only the leading front is drawn. For the sake of comparison, the corresponding self-similar curves, denoted by “ss”, are plotted too. The double arrow and Of represent the local Debye length and the ion velocity, respectively, at the position of maximum electric field.

376

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A liydrodynamic approach

In the following we shall prove that an intermediate stage exists indeed in which quasi-neutrality holds, but in which the self-similar behaviour has not yet been achieved. This stage will be characterized by long wave-length oscillations. The starting point for our analysis is the scalar wave-equation in the quasi-neutral limit given by V+a~(hIV)=0,

(10.3)

which has been derived in section 6.1 and assumes ‘y = 1. V is the specific volume (V= 1/n), and ~j the Lagrangian mass variable defined by eq. (6.1). The self-similar solution of (10.3) is given by

VJ~’vj,t) = _t/*q —1 ~‘q/t~0,

(10.4)

(see section 6.1, eq. (6.16)). We now look for a perturbative solution of (10.3) around V05 making the ansatz: V=V~5(1+w), where wi

(10.5)

1. Inserting (10.5) into (10.3) we get after linearization 2a a~w) 2 a1(t 1w) = a~(~ By means of the transformation =~

.

u=—~7tw

(10.6)

(10.7)

(10.6) becomes t2u=~2u”,

(10.8)

where dot represents a 1 and prime a~,respectively. Equation (10.8) can be solved by separation of variables u(ij,t)=f(t)g(~).

(10.9)

The symmetry in (10.8) implies that f(t) and g(~)are governed by the same differential equation 2yldx2=(A1x2)y, (10.10) d which is called the Euler differential equation (y(x) = f(t) or g(~),respectively); A is the separation constant. For 4A + 1 <0 the solution reads [82]: f(t) =

‘s/j

[a 1 cos(s ln

t) +

a2 sin(s ln t)]

g(ij) = -s/~i~[b1 cos(s ln(—q))

+ b2

,

sin(s ln(—~))],

(10.lla) (10.llb)

C/s. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

377

where a~,b~,i = 1,2 are constants, and s is given by s=~Vi4A+1i. To simplify the discussion, we set a2 V=

+

= b2 =

0. With this, V in (10.5) becomes

cos(s ln t) cos(s ln(_~))],

(10.12)

where c is constant. For fixed t the oscillatory pattern in ~j is evident, a fact which transforms to the other quantities and to the Eulerian space as well. Equation (10.12) is valid for large values of t and Intl. Excluding the front given by ~ 0, we see that the oscillation damps out in the course of time. Time-asymptotically the self-similar state is reached. To our knowledge, eq. (10.12) represents the first consistent, analytic expression in the expansion problem which goes beyond the self-similar treatment. Since charge separation is neglected it applies to the gas dynamical expansion as well. The consequences of this observation can be stated as follows: 1. There exists a state of intermediate asymptotics [54] in which the solution no longer depends on the details of the initial and/or boundary conditions but is far from being in the limiting state. 2. This limiting state is self-similar and is stable. 3. The final approach to the self-similar state is not affected by charge-separation in the plasma case.

11. Comparison with experiments Despite the simplicity of the considered ion expansion model comparisons with laboratory experiments are possible. Plasma expansion experiments can roughly be divided into two not necessarily’~omplementary classes, namely into i) experiments dealing with vacuum arcs, exploding wires etc. in which a metal plasma produced at the cathode surface expands towards the anode and into ii) experiments that are more or less designed for the study of the quasi one-dimensional expansion of an initially step-like plasma (“diaphragm problem”, “dam-break problem”). Ion acceleration experiments in laser produced plasma are purposely omitted here. Their explanations have to take into account among others the irradiation field, suprathermal electrons, several ion species as well as possible nonplanar geometries, effects that have been ignored in our model. A review on the first set of experiments has been recently given by Jüttner [83]. It is generally believed that a major mechanism for the acceleration of the ions in vacuum arcs is the electron pressure gradient [21,22,24]. The ions easily gain energies exceeding the applied voltage, a process which appears to be enhanced by the Joule heating of the arc current or by anomalous heating processes [84]. Cathodic erosion creates craters at the cathode surface, the cathode spots, giving rise to a continuous quasi steady-state process. (Mathematically one encounters here a somewhat different boundary value problem. Instead of a rarefaction wave propagating towards the unperturbed plasma described asymptotically by the time-dependent self-similar solution, there is a quasi-stationary solution being due to a plasma source at the cathode surface i.e. at a fixed boundary.) As pointed out by Mesyats [85],

378

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

however, it cannot be ruled out that the high velocity ions with U (1—2) X i0~rn/s are ejected by explosive surface processes, resulting from a process which seems to be more relevant for exploding wires. The second set of experiments in which this later non-plasma dynamical process is absent allows a closer contact. Korn et al. [26] investigated the evolution of a large density discontinuity in a single-ended 0 machine and observed essentially the self-similar motion corrected by a finite initial ion drift velocity U(x, 0) and by a finite ion temperature, T 1 = Te. When the electron temperature was raised by -=

microwave heating at cyclotron resonance, the spatial profile was observed to broaden, as expected. A peak of fast ions was, however, not mentioned. Our numerical results indicate that the finite ion drift velocity, a characteristic quantity in his experiment, suppresses particle bunching (see section 9 and fig. 13). A group of fast ions has been recorded by Tyulina [25]who measured the transverse expansion of a pulsed, nonmagnetized arc. In this experiment the existence of a group of ions preceding the plasma was seen with energies on the order of ten times the estimated 2, ~ electron being thethermal initial energy, thermal corresponding energy of the to maximum ion velocities of about 3c5, c~ = (kTeo/mi)” electrons. Similar maximum ion velocities were found in [22—24]. A series of detailed measurements have been presented more recently by Eselevich and Fainshtain [30] to which we draw the attention now. In this experiment the plasma was produced in a prechamber at the entrance of a big vacuum chamber by ionizing pulse-fed gas. During plasma expansion the density and electron temperature were nearly constant in this section of the chamber. In spite of the modified boundary conditions, the results bear resemblance to ours in several respects, as follows. The most remarkable one was the measurement of fast ion velocities of 15 c 0 at large distances from the source, corresponding to an energy E 100 kTeo. The acceleration took place under isothermal conditions initially, later on Te decreased in the front region. These high ion velocities are contrasted with the relatively low ion front velocities in the earlier experiments [22—25].One interpretation of the latter is based on the small dimensions of the used vacuum chambers, which were not larger than 0.2—0.4 m. Consequently, in these experiments only a limited time interval was available for the expansion; the fastest ions could only be accelerated to velocities of about 3c~(see, e.g., the plateau regime in fig. 18, the curve for y = 1). In the experiment of Eselevich and Fainshtain, however, the vacuum chamber had a length of 2.0 m. By this it was possible to follow the expansion of the plasma over larger distances and, thus, also over longer times. Therefore, the time evolution of the front velocity could reach the third phase in fig. 18, the curve for y = 1, in which the velocity monotonically increases and in which it assumes values much larger than 3c5. Another reason for the reduced ion front velocities in the experiments [22—25]may be attributed to the fact that the electrons do not behave isothermally during the expansion. In this case the ion front velocity time-asymptotically approaches the values given by eq. (10.1) which are in accordance with our numerical results (see fig. 18, the curves for y ~ 1). More than a qualitative explanation of the expansion process cannot be expected within the framework of our simplified model. As an example, the plateau regime found in the numerical calculations (see section 10.2, fig. 18) has not yet been detected experimentally. Furthermore, detailed experimental data concerning the thermodynamic behaviour of the electrons in the initial state near the generation point of the plasma and during the expansion are needed, to make estimates for the available ion front velocities. Eselevich and Fainshtain [30] compare their results with well-known kinetic [33,46] and hy-

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

379

drodynamic [42] models for the collisionless expansion of a plasma into vacuum. Especially, the hydrodynamic model appears to be quite suitable for the description of the dynamics of the accelerated ions, because the measured effective ion temperature rapidly decreases in the direction of the expansion. Consequently, the thermal energy of the ions can be neglected compared to their kinetic energy (see also [33]). The comparison of the collision times and of the mean-free paths with the characteristic time- and space-scales of the experiment shows that the effect of binary collisions and of collisons with multiple small angle deviations on the expansion is negligible within the temperature and density regions given by Eselevich and Fainshtein: 3. (11.1) 1-4eV, nei = 107_lOb cm For collisions of electrons and of ions with neutral particles Eselevich and Fainshtain [30] (1980) give a mean-free path of more than 1.0 m. The collision times Te and r~as well as the mean-free paths for electron—electron and ion—ion collisions are obtained from (see [86]): Te

=

=

2-10eV,

T1

=

3.44 X 10~T~2/n~A, 7

3/2

2

le

= UtheTe

(11.2)

1/2

2.09 x 10 (T~ /Z n 1A)p.and ,Uthi = 9.79 ii = UIh1T , (11.3) 7V~cm/s x i0~\/~7~ cm/s denote the thermal velocities of the where Uthe = 4.19 x iO electrons and ions, respectively. Z, ~ and A are the charge number, the fraction of ion to proton mass, =

=

m

1Im~,, and the Coulomb logarithm, A 10—12, respectively. The temperatures Te and T1 are expressed in 9eV.cm’3, For the single ionized argon gas (Z = 1, ~ = 40) used in the experiment one gets with, Te = 6eV, T e.g., ~ei = i0 1 = 1 eV and A = 15 the following values: Te

sxs4 X 1O”~s, 2s, 1.3 x 10

Uthe Uthi

10 cm/s, 1.5 X i05 cm/s

(11.4)

and le~4x104cm,

l

3cm. 1~=2X10 The calculated collision times Tei in (11.4) are much larger than the flight times of the particles through the vacuum chamber being of the order of (3—5) x i0~sfor the ions. The mean-free paths lei exceed by far the maximum flight distance of 2 x 102 cm. For the time evolution of the ion front velocity Eselevich and Fainshtain find a logarithmic dependence being in good agreement with the numerical results of the isothermal expansion model of Crow et al. [42],at least with respect to the slope. However, as can be checked by means of ref. [30]fig. 4 (1981), the numerical value of the slope does not correspond to the given scaling law, ref. [30] eq. (11) (1981) (see also eq. (3.35) in this review). Measurements of the spatial structure of the ion front in the experiment of Eselevich and Fainshtain have yielded two different expansion regimes in which the front width 4 can be either larger or smaller than 10 AD! [30]: AD! is the local Debye length at the front. For 4 ~ 10 AD! Eselevich and Fainshtain state

380

Cli. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

a quasi-neutral behaviour of the front, where the density falls off almost monotonically. In the other regime, where the condition 4 a 10 AD! is satisfied only at the plasma source, the density steepens during the expansion, and the front width 4 varies in the interval AD! ~ 4< 10 AD!; furthermore, a density hump appears in the ion density profile (see ref. [30] fig. 1 (1979) and ref. [30] fig. 5 (1980)). This behaviour is typical for the presence of charge separation effects at the expansion front. Our numerical results in section 5 and in section 9 offer a possibility to explain the front structures found in the experiment of ref. [30]. In section 5 it has been stated, that the steepening of the density profile and the subsequent bunching are driven by a spatially decreasing velocity profile resulting from the hump-like electric field of the initial state (see e.g. fig. 2). Steepening and bunching are particularly pronounced, if the initial velocity of the plasma is low (see section 9 and figs. 13, 14), and thus the ambipolar field dominates in the expansion process. Such conditions could have been present in the experiment of Eselevich and Fainshtain in those cases where the decrease of the front width from initially 4 ~ 10 AD! to AD! ~ < 10 AD! could be observed. In the other cases in which always 4 ~ 10 ADI was measured during the expansion it can be concluded that the electric field at the plasma source and / or the initial velocity of the plasma was spatially increasing and therefore, no significant steepening could take place. However, similar to the interpretation of the behaviour of the maximum ion velocity at the beginning of this section, we are left to a qualitative description of the occurring processes. For a unique clarification of the space-time behaviour of the front observed in the experiment of Eselevich and Fainshtain detailed information about the plasma at its generation point are missing. This holds true for any other experiment, too. ~

12. Summary and conclusions In this review we tried to summarize some aspects of the hydrodynamic evolution of a plasma expanding into vacuum. A primitive form of the description was chosen neglecting important effects such as the finite size of a plasma, magnetic fields, realistic collisions, two-electron temperatures, several ion species, and so forth. On the other hand, the chosen model turns out to be by no means trivial. It allows the study of the plasma expansion as a well-defined initial and boundary value problem. The study includes the complete time evolution, starting from the very beginning where the plasma distribution is step-like up to the final self-similar state. A variety of dynamical structures have been found and explained. In the dissipationless case the occurrence of wave breaking in connection with ion bunching is one of the most striking features. An analytical confirmation of this phenomenon was given by deriving and solving a scalar wave equation. The latter describes the whole ion dynamics inclusively charge separation. The observation of wave breaking was rendered possible by using a dissipationless Lagrangian scheme with moving boundaries. In addition, an effective Poisson solver was developed involving the smoothest dependent quantity 4’ only. Long term calculations could be performed by introducing a Navier—Stokes viscosity. It stabilizes ion wave collapse. As a relic of the nonlinear steepening process, a bunch of fast ion appears. In general, three dynamical phases have been identified. The first one is the initial acceleration phase. The second one is characterized by the presence of the stabilized ion bump, the velocity of which is approximately constant (“plateau regime”) and depends on y, the electron polytropic exponent. The third one is the final acceleration phase starting with the debunching of the ion bump. It was shown that the whole

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

381

plasma, including the front behaviour, tends to the self-similar state. In addition, a state of intermediate asymptotics has been detected. Although some of the features seem to be in accordance with experiments, refinements are needed to cope with realistic expansion processes. Ionization, recombination, Coulomb collisions, neutral particles, electron heat conduction, and finite energy reservoir are some of them. There have been efforts to incorporate some of these effects in plasma expansion models [28, 29, 40, 41, 51, 52, 87, 88] inclusively the investigations dealing with the ion blowoff from laser-produced plasma [6—11,27—29, 45—52, 89—92]. A reasonable understanding of realistic plasma expansion is, however, still waiting for. One reason is the intricate nature of the involved intrinsic processes, another one the occurrence of nonlinear dynamical structures, some of which having been elucidated in this report.

Acknowledgements The authors would like to thank Prof. Dr. K. Elsãsser for fruitful discussions on the topic of this review. Furthermore, the authors thank Mrs. H. Schamel for translating the text and Mrs. B. Malik for typing the manuscript.

Appendix A: Bunching and wave breaking in nonlinear Langmuir oscillations The dynamics of nonlinear Langmuir oscillations can be described within the framework of a macroscopic model. The hydrodynamic equations for the electrons are completed by Poisson’s equation. The ions are assumed to form a fixed neutralizing background of uniform constant density. Furthermore, the electrons are considered to be cold, so that the pressure is neglected in the electron momentum equation. From these assumptions the following set of equations is obtained (Davidson [62]):

atne + ax(neUe) atUe +

Ue

=

axUe

=

4rre(ne

=

0

—(eIm~)E —

n0).

(Al) (A2) (A3)

Equivalent to Poisson’s equation (A3) is the Ampere’s law in the electrostatic approximation —

41ren1,t.le = 0.

(A4)

Within the context of eqs. (Al), (A2) and (A4) Poisson’s ‘equation may be considered as an initial condition. By virtue of the continuity equation (Al) and eq. (A4), eq. (A3) remains true for all time if true initially. The linear analysis eqs. (A1)—(A2) Assuming shows thatathe only modes of oscillation arefor at the the electron 0pe of = (41Tnoe2Ime)~”2. sinusoidal initial perturbation electron plasma frequency ~ density, one obtains solutions corresponding to standing waves (wave length 21T/k, frequency Wpe)’

382

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

which retain their sinusoidal wave form. In the case of large wave amplitudes, however, where the linearization procedure is inapplicable, strong distortions of such wave forms have to be expected. This is due to the dependence of the solution on higher harmonics in kx. Indeed the concept of distortionless oscillations is completely untenable in the investigation of eqs. (Al)—(A4), if large-amplitude perturbations are admitted. An exact solution of eqs. (Al)—(A4) is rendered possible by the introduction of Lagrangian variables [62] corresponding to the ion dynamical case (see section 4.1). The set of equations transformed in this manner reads: (AS) (A6) fle(~,T) ~X(~,T)=fle(~,0)

(A7)

with mast and ~‘=X~j~dT’Ve(~,T’). Equations (A3) and (A4) can be combined to give E(~,T)

=

4~fl0~U~(~ T).

(A8)

Differentiating eq. (A6) with respect to m and using eq. (A8) one obtains for oscillator equation

Ve( ~,

m) the harmonic

(A9)

~Ue(~,T)+W;eVe(~,T)=0.

The advantage resulting from the transformation of the set of equations (Al)—(A4) is evident. Equation (A9) is a linear scalar wave equation, which can be easily solved compared to the original set in the Eulerian picture. The Lagrangian position variable 4 appears only as a parameter. From the general solution of eq. (A9) Ue(~, r) = V(~)

cos

T

X(~)sin Wpe r

+ ~

(AlO)

the remaining quantities can be derived by integration or differentiation of eq. (AlO): x(~, m) =

~+

sin

+ 0pe

E(~,r) ~e(~’



‘ “

=

r) =

Wpe

~

[V(~) sin

0) [1 +

X(~)(1

m



t0pe



cos Wpe m)

X(~)cos

(All) (A12)

Wpe T]

~

V’(~)sin

Wpe T +

represents the differentiation with respect to

X’(~) (1 —cos ~.

Wpe T)]

;

(Al3)

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

383

The functional dependence of V and X on ~ in eqs. (AlO)—(A13) is related to the initial conditions for the velocity and for the electric field as follows: V(~,O)=Ue(~,O),

X(~)=-

2

(Al4)

E(~,0).

mewpe

In addition, X(~)is related to the initial density Tte( ~, 0) through Poisson’s equation at m X’(~)= fle(~,0)~’Tbo —1.

=

0, i.e. (A15)

For the transformation from Lagrangian to Eulerian variables, i.e., the determination of ~ as a function of x and t from eq. (All), the initial conditions V( ~) and X( ~) have to be specified explicitly. Generally, this transformation entails the solution of a transcendental equation which can be established by numerical methods. The condition, that the electron density should be positive and finite for all time, restricts the class of initial-value problems which may be treated by the procedure outlined above. Additionally, the following inequalities have to be fulfilled: ne( ~, 0)> n0/2 1 Wpe

/

(A16) (•~

V’(~)i<(2 ~ flo S’ \

0~ —1)/

\!/2

(A17)

.

Mathematically, if inequalities (Al6) and (A17) are violated for some ranges of ~, the transformation from the Lagrangian to the Eulerian variables loses its uniqueness (see eq. (All)). Physically, the violation of (A16) and (Al7) leads to the development of multi-stream flow. In particular, inequalities (Al6) and (Al7) exclude electron trapping within the context of the present cold-plasma model [62].To explain this statement, consider the initial conditions X(~)=O,

fle(~,0)fo,

V(~)=U0sink~.

(A18)

From eq. (A 10) and inequality (Al7) it follows that the maximum electron flow velocity is restricted during the evolution of the system, i.e.

k’eimax =

VOl

<

iWpe/ki

(A19)

.

Since the plasma is cold, no electrons are moving exactly at the phase velocity Iw~~/ki of the fundamental wave, i.e., there exists no resonance between the electrons and the fundamental wave, and therefore trapping does not occur. Next we present an example taken from ref. [62], which demonstrates the essential processes connected with the transformation to Eulerian variables. The initial conditions are specified as follows: ne(~,0)no(1+4c05~),

4~<~

(A20)

and Ue(~,O)=V(~)=O.

(A21)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

384

Inequalities (Al6) and (Al7) are satisfied for this case, and one obtains from eqs. (A10)—(Al3): Ue( ~, T) =

4 sin k~sin

—~

m E(~,m)

=



m

(A22)

2

e

sin k~cos

pe

m

(A23)

I + 4 cos k~ ~

(A24)

with kx

= k~ +

a(T) sin k~,

a(r)

=

24 sin2(wper/2).

(A25)

In (A20) the condition ~I< ~ guarantees that a(r)i < 1 and that the inversion of (A25), i.e. x(~,m)—~~(x, t), is unique. For small amplitudes, i.e. I~I ~ 1, it is clear that the Eulerian and the Lagrangian variables are approximately the same, x 4. In this case the results correspond to those which follow from the linearization of the equations of motion (Al)—(A4). For larger values of I~I,however, the solution has to be determined numerically from (A25) at different times. Then Ue~ E and ~e’ in Eulerian variables, are deduced from eqs. (A22)—(A24). Solutions of this kind are presented by Davidson [62]for 4 = 0.45 (see ref. [62] figs. 3.1 and 3.2). The xt-dependence of the solution is periodic with periods 2ir/k and 2~m/w~. The most interesting property of the solution is that the electrostatic forces cause a strong local increase of the electron density in the vicinity of —

x(2n+l)ir/k,

n0,±1,±2,.

.

(A26)

.

This phenomenon is called “bunching” (compare [62]). The density maximum,

lie = 5.5

n

0 is reached at t = ir/w~.Later on the system returns to its initial state. The corresponding electric field exhibits a steepening of the initial wave sin kx with any increase in the maximum amplitude. The electron velocity, however, does not steepen; for t = (IT/Wpe) Ue is even equal to zero. The solutions presented by Davidson [58] have a certain similarity to the solutions for the ionic system shown in fig. 7. However, in the case of ion dynamics, the velocity steepens, too (see fig. 7). Such a behaviour can also be obtained for nonlinear Langmuir oscillations, if one uses a more general formulation of the electron motion. This formulation, allowing the separation into a boundary- and an initial value problem, has been established by Kalman [73]. Different to Davidson [62], Kalman did not use the initial coordinate ~ as a Lagrangian variable, but the so-called Lagrangian mass variable n(~I’in Kalman’s work). However, there is a unique relation between ~ or x and n:

:=

Jn5~’,

0) d~’

fle(X’, t)

dx’.

(A27)

The transformation (A27) automatically satisfies the continuity equation

an/axl, so that atne

=

fle(x, t),

+ ax(neUe) =

an/ati~= 0 is valid.

~eUe~

(A28)

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

Additionally, from (A28) it follows that the total time derivative of dn/dtO,

385

n vanishes: (A29)

i.e. n = const. along the stream lines. Therefore, n is sometimes called stream function. Inserting the transformation (A27) and eq. (A29) into the set of equations (A1)—(A4), one finally gets again the scalar wave equation for the velocity Ue: (A30)

2Ue/at2+w~eUe=o.

a

Different to the wave equation (A9) Ve now depends on n and t and not on ~ and r(=t). In Kalman’s work [73], however, a scalar wave equation for the velocity is not deduced but a corresponding equation for the space coordinate x, which contains a free function of time. The solution of this equation yields the velocity, the density, and the electrostatic fields in dependence on n and t. For an explicit representation of the solution the boundary and initial values have to be specified. Defining, e.g., a boundary value problem of the kind Ue

= Ue

0(l + ~.t

.1

sin tot)

for x = 0

= ~eü~’eü

(A3l) (A32) (A33)

E=0

where ~ is the relative depth of the velocity modulation and to is its frequency, one obtains solutions of the problem, the structure of which is governed by the parameter Co/COpe.

(A34)

For a = 2, e.g., the bunching of the electrons is unlimited (the electron density becomes singular!), and the electric field and the velocity steepen indefinitely with a reversed sign of their gradients (see ref. [73] fig. 1). This property of the solution in the electronic case corresponds exactly to the properties of our solution shown in fig. 7. It is, therefore, obvious to investigate the nonlinear ion motion with a formalism analogous to that used by Davidson [62]and Kalman [73]. The scalar wave equation for the ionic system deduced in section 6.1 is, however, partial and nonlinear. Therefore, it can only be treated by perturbative means (see section 6.2).

Appendix B: Spatial discretization of Poisson’s equation According to Potter [66] a double mesh point net is used for the spatial discretization of the hydrodynamic ion equations in the Lagrangian picture, eqs. (4.9)—(4.ll), respectively (8.1), and of Poisson’s equation (2.24). The cell boundaries x1 and the velocities U. are defined pointwise on a net J, 1 ~j ~ J + 1; J denotes the total number of cells in the finite integration interval having the length 2L, —L ~ x ~ +L, at t = 0. The properties of the density and of the electrostatic potential, respectively field are distributed within the cell. Therefore, the mean points of each cell are defined on a second net I,

386

1

Ch. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

i ~ I(=J); each cell i has the width

4,

=

x1÷1 x. —

(Bi)

.

The relation between the f-net and the I-net is given by x~=

~(x1~1 x1).

(B2)

+

~.x

At t = 0 the grid points have equal distance: 4, = = const. For our purpose the choice f = 500, L = 200, and thus, ~x = 2L/J = 0.8 has been sufficient (see also [56]). Since Poisson’s equation is only spatially discretized, the time variable is not explicitly considered in the following calculations. The discretized version of the linearized Poisson’s equation on the I-net reads —

(dfl~/d4’1~)~4’~ = ~(dfle/d4’o)i 4’~ + n~(4’~) n,.

(B3)



4’ in

For the second derivative of =

4’~,

-

k=

(B3) at the grid point i one obtains (B4)

i +

xkxkl

with 4’j+1 —



i+l

-

,



x~





4’~-~-

(



and xk

=

~(x 1~1 + x~),

Xkl =

(B6)

~(x~ + x~1) .

Inserting (B5) and (B6) into (B4) the resulting equation yields a three-point formula, i.e. 4~ is completely determined through the discrete values of 4’ at the grid points x, x, and x. ~ By this one gets for (B3) a set of linear equations with a tn-diagonal coefficient matrix which can be formulated as follows ~,

~

2~i~I—l

(B7)

with ~

= —

1

a,,

= —L

= /

2 x~_1j~x~ x~1)

(B7a)



2 /dn\1 + ~—~) ] —x,)(x, —x~1) dçb0

(x,÷1

2 \1 \ — x1_1,~x~÷1 — x~)

(B7b) (B7c)

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

387

in (B7) o is identical with the right-hand side of (B3). The quantities fle(4’01) and (dne/dq5o)i are given by (4.14) and (4.15), where 4’~denotes the values of 4’, from the previous iteration step (see section 4.2). The representation of (B3) by (B7) cannot be transferred to the first and to the last row of the coefficient matrix because in both rows only two instead of three matrix elements are available. Note that a second derivative has to be resolved by three grid points. The missing matrix elements are defined by the boundary conditions being deduced in section 4.3, eqs. (4.37) and (4.42). These equations are differential equations of first order and are approximated in the discrete case by two-point formulae. The first row of the coefficient matrix (i = 1) contains the boundary conditions at x = — L (see eq. (4.37)): (B8)

all4’l+a124’2=o-!=O

with a11

—4’~

=

a12 =

4’o2 +

ô

(B8a) — +

1

whereas the last row (i = +

x3—x2

(4’~~ —

a114’1 =

:=



is filled by the boundary conditions at x =

I)

-

,

+L

(B8b)

(see eq. (4.42)): (B9)

0~~

with a11..1

=

—[1 —f11...1(x1

a11

1

+

=

=

f11



(B9b) g1



~

(x1 — x1...1)

(B9c)

1~~’2 + [n {[~e(4’oi)]

1~2}

(B9d)

~

{[n~(4’o1)]~

0(4’01.1)]

~

2+

g 1,1~1:=

(B9a)

1(x1 — x1..1)

[f1~...1(4’~1 + ~

fj,1-1

x1.1)]

[ne(4’oji)]~2.

(B9e)

The numerical solution of the set of linear equations represented by eqs. (B7)—(B9) is accomplished through a Gaussian elimination procedure adapted to the tn-diagonal band structure of the coefficient matrix (see [93]). Since generally Poisson’s equation (2.24) is nonlinear (exception y = 2) and, therefore, can only be solved iteratively, the outlined procedure has to be repeatedly used for each time step in order to determine the new values of the potential. The maximum relative error defined by (BlO)

Ch. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

388

is used as a convergence criterion for the iteration; p numerates the iteration step. In the case of convergence 6.,~has to decrease monotonically at least from a certain iteration step on, so that ~6 can be achieved. Appendix C: Implicit difference scheme for the equations of viscid ion motion For the time discretization of the equations of ion motion (4.9), (4.10) and (8.1), the time t is divided into finite equidistant increments ~t, i.e., t = n~t=: t~,n = 0, 1,2 In order to guarantee an accuracy of the order O(~t)2it is necessary to define each cell boundary at times n, whereas the velocities U. of the cell boundaries are defined at the intermediate time steps n ~ [66]. By this one obtains in the Lagrangian picture the following difference scheme on the double grid point net (see appendix B): —

U~’2) /~t = ~(E~



+

E~

2 +

(1



a)(U”)~”2]

(Cl)

1)+ v[a(U”)~’

(C2)

XX~+~tU~~ n+1

n~

n

n

n+1

fl

n~(x

n+1

),

1+~—x1)/(x1+1 —x1

(C3)

where E’ in (Cl) follows from !

fl

~f

fl

‘~

fl

—~I(4’~+~ — 4’1)/(x~+~ —x1)+ (4’,-

~

fl

‘~

4’~1)/(x~ —x11)]; (C4) 4’~results from the iterative solution of Poisson’s equation (see appendix B), and x7 is given by (B2). Except the calculations for fig. 7, being performed with z~t= 0.005, z~thas been set equal to 0.025. The spatially centred formulation of the second derivative of U in (Cl) is given by E1

=

(v

)j

((U

m )J+1-2

(v



,

m )~~1-2)



m rn /(x1+~12 — x1_117),

m = n±~,

(CS)

where m

x1~12 =

~(x1~1 x1) 1

~n

+

~

(C6)

and in

m

in

(U )J+1/2 = (U1÷1— U1

m

in

)/(x1~1 x1 —

),

in

x1

=

1

fl

1(x1

n ±1

+ x1

).

(C7)

In (Cl) the quantity a denotes the implicitness parameter, 0 ~ a ~ 1. For a = 0 the difference is 2.Inscheme the case explicit, because the right-hand side of (Cl) only contains U of the previous time step t”’~’ a ~ 0 U of the actual time step t~/2 comes into play, and the scheme is implicit. Special cases are a = (half-implicit) and a = 1 (full-implicit) (see section 8). If the viscosity constant i.’ is equal to zero, the dissipationless explicit Lagrangian scheme is obtained which has been used to get the results presented in section 5. Similar to the discretization of Poisson’s equation in appendix B the spatial discretization of (Cl)

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

389

leads to a set of linear equations with a tn-diagonal coefficient matrix: 2, b~1U~1+ b~V~ + b~~1U~~1 = r’~

1

n

+

2~j~ J.

~,

(C8)

The coefficients in (C8) and the right-hand side u~”~112 are obtained by inserting (C5)—(C7) into (Cl): b~ 11= —s/4~1j1 b~1= 1 ~

+

~

a

(C8a) (C8b)

~

~‘~J+I,J!

=

s = 2~At

~

2=

(E~ +

(C8c)

~

2+

E~

r’~”

p

At (1— a) (U”)~2,

(C8d)

1)+ v~ =

x~ 1—

=

,

x~+~ — x~,

4~-~ =

— x~1

.

(C8e)

The first and the last row of the coefficient matrix are specified by the boundary conditions in the finite integration interval (see section 4.3). At the left-hand boundary, x = —L(j = 1), the discretized version of eq. (4.33) yields: (C9)

b’11V~+ b~2t4= r~’; b~1 1 +2Ata-’,h/(x’~—x~)

(C9a)

b~2= 1 —2 At a\/5~/(x’~ — x~) 1+ U~~1 +2At(l — a)-~~/y (Ut’ ~ = V~

(C9b) (C9c)



with I from (C8). At the right-hand boundary, x eq. (4.44)) is used:

v~’)I(x~~ —xv’) =

+L(j =

J +

1), the spatial mean value of t~=

E

(see

+ t)~,)= E

(ClO)

1.

From (ClO) it follows: I

/

~

I / n,n—1/2 + b~+1~+1U1+1 = r~~1

,

(Cll)

where b1~1~ = 1

(Clla)

= b1~11~,

and 2= V~1+ ~ r~11U

+

2 AtE~ .

(Cllb)

C/s. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

390

In the implicit case, a 0, an iterative procedure is required for the determination of the velocity 2 in (Cl) by means of (C8), (C9) and (Cll), which can be understood as follows. Assume that f, U~U U”112, E~and n’ are known at the time step n. In order to calculate U’’~2 from (Cl), one needs, according to (CS )—(C7) x’ however, x’ + is not known, since for its determination the knowledge of UnL(2 is necessary (see (C2)). Consequently, one has to iterate between (Cl) and (C2) at each time step. As initial value for iteration the values of x from the previous time step are used. Typically two or three iterations are necessary to reduce the maximum relative error 6~below i07 6~is defined by 1;

:=

maxi(U~



(C12)

U~)/(U~’ + V~)I,

where p denotes the iteration step. An expression similar to (C12) holds for xc. At each iteration step the set of linear equations consisting of (C8), (C9) and (Cil) is solved by means of the Gaussian elimination procedure mentioned in appendix B.

Appendix D: Conventional stability analysis of the discretized set of equations In section 8 the assumptions for the derivation of a conventional stability condition of the implicit difference scheme have been discussed, and its accordance with the stability behaviour of the numerical solution has been investigated. A crucial point for the derivation of the stability condition is the use of the pseudo-sound velocity (see section 3.2 and eq. (8.5)) and the assumption of an appropriate linearization of the ion equations (see e.g. eqs. (8.3) and (8.4)). Thus, the equations (4.9), (4.10) and (8.3) are discretized with respect to the Lagrangian coordinate ~ to give: (Dl)

(x~—x~)!At=U~2 V 0(x~~1 —xc)/h=Vc+112 n+1/2

(D2)

n—1/2

U-

C

U-

At

= ~

(V~÷3~2 —

V~1~2 + V7±112 2 V~~’~2



v’~--’

+~

[a ~ I

+BLa

h

Vc_312)

+(l—a)

‘‘

n+1/2 V-~~ —2V-n+1/2 +V-n+l/2 1 ‘2

v”~2 v~112 j+1

+(l—a)

j-I ] n—1/2

U-~1

n—!/2

—2U-

n—1/2 +U-1

1

~2

(D3) where V0 = 1 /n0(fl denotes the specific volume at t = 0, and h(=A~)is the equidistant space-step; a is the implicitness parameter (compare (8.2), resp. (Cl)). The quantities A, B, C in (D3) are defined by (8.4) and are held locally constant, as well 1/~in (D2). By means of the discrete time differentiation of (D2) and by using (Dl) the space variable x can be eliminated: V0

At -~-

n—! /2

n—1/2

(v~÷~ — U1

n ) = VJ÷1/2

n—i



VJ+!/2.

(D4)

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

391

(D3) and (D4) represent a closed set of linearized equations, the stability of which is investigated in the following. Inserting the Fourier transformation

f~=

j~exp(ik~

1)

(D5)

into (D3) and (D4), where obtains:

denotes the discrete values of ~ and f stands for U and V respectively, one

(w

= y 12;1 +

+

ry)

(D6)

~n-1/2

(D7)

2. rC~

+

= ~n-1

In (D6) V~has been expressed by (D7). For the coefficients r, y, u and w, resp., one gets: r = 2i I’~~ sin(~)

(D8a)

y

(D8b)

= iC

~ cos(~) sin kh

u=l+2(l—coskh) C 2—2i sin khC1 w

=

1 —2(1



cos kh)

~2 +

2i sin kh

(D8c)

ë~

(D8d)

with /C1\

/

1=1

I

a

~AAt I—

\l—a/

and

2h

/

a

\BAt ~C2~ \l—a/ h

/C2\ I



I=I

The eigenvalues of the coefficient matrix of the linearized set (D6, D7), i.e. (~

(D9)

(w+ry)Iu)

yield the stability condition. The characteristic equation for the eigenvalues of (D9) reads 2— A(u + ~

ry) +

=0.

(Dl0)

A

In general its complex solutions are 2

A! 2

=



4uw)

(Dll)

(L ±VL

with L=u+w+ry.

(D12)

392

Cli. Sack and H. Sc/same!, Plasma expansion into vacuum — A hydrodynamic approach

The system is stable, if A11

~ 1, j = 1, 2

(compare [65, 66]). L and the square root in (Dll) can be

2a) + iT(1



2a)]

+

Y~ T2

formulated as follows L

=

2[1

2 Y



Q



1(l



(D13) (D14)

2—4uw=2\[X~iY \/L

with X= Y=



—2T

2Q2 + 2Q2

~“1 (1

[Q2(l

+



2a)



2a)

(Dl5a)



Y 1]

(DlSb) (DiSc)

2 = ~p2 sin2 kh

Q

BAt (l—coskh) T= BAt af(ln p2

=

(D1Sd)

sin kh

~)

(DlSe)

c2 (1/a At/V/i)2,

(DlSf)

where c2 is given by (8.5). From these quantities one gets for A! 2 in (Dll):

[L ± \/~(Vi~+ X + iV7~ X)] u~,

= ~

(Dl6)



2 >0; U is given by (D8c), and u’~is the complex conjugate of u. From A where R = + Y 1 one deduces after a lengthy calculation the relation

2 ~1

2I

(D17)

0~R2~2Y~—X.

After a further quadrature it follows O~Y2 ~4Y2(Y~



X).

(Dl8)

(Dl8) is the condition for numerical stability we have looked for. According to the sign of c2 two cases have to be distinguished. For c2 > 0 one obtains from (Dl8) and (DlSa):

Q2 [T2 (1 Q,



2a)2 +

Y~] [2Y~ ~



2Y

2 1 (1



+

Y~)].

(Dl9)

2a) (T

Y

1 and Tare defined by eqs. (l5c—e). From this the cases a = ~, inequality (8.6), and a =0,1, inequality (8.7), are deduced, where in the latter case cos kh = 0 (kh = ir/2) has been used.

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynamic approach

393

In regions of increasing density, i.e. c2 <0, (D18) reads:

1Q21 [T2 (1

2 + Y~)]. (D20) 1 (1— 2a) (T In comparison to (D19), (D20) has a negative sign on its right-hand side. This means that the set of discretized equations is numerically unconditionally unstable, if the bracket on the right-hand side of (D20) is positive. Since ~! is always greater than zero, this is the case, if a is set equal to one, i.e., in the full-implicit scheme. However, such an instability has not been detected in our calculations. —

2a)2 + Y~]~ —[2Y~ — 2Y

References [1] U. Samir and N.H. Stone, Bericht der “XVI International Conference on Phenomena in Ionized Gases (ICPIG)” Düsseldorf, FRG, 1983; p. 84; K.H. Wright Jr., N.H. Stone and U. Samir, J. Plasma Phys. 33(1985) 71. [2] i.E. Borovsky, MB. Pongratz, R.A. Roussel-Dupré and T.-H. Tan, Astrophys. J. 280 (1984) 802; T.-H. Tan and J.E. Borovsky, Preprint LA-UR-84-754, Los Alamos National Laboratory, Los Alamos, New Mexico, USA, 1984. [3] R. Behrisch, in: Proc. XIth Intern. Symp. on Discharges and Electrical Insulation in Vacuum, Vol. 1, Ostberlin, GDR, 1984 p. 123; L.M. Baskin, AL. Gromov and G.N. Fursey, ibid., p. 165. [4] PB. Parks and R.J. Turnbull, Phys. Fluids 21(1978)1735. [5] CT. Chang, LW. Jorgensen, P. Nielsen and L.L. Lengyel, NucI. Fusion 20 (1980) 859. [6] K.A. Bruckner and S. Jorna, Rev. Mod. Phys. 46 (1974) 325. [7] E.J. Valeo and Ira B. Bernsteiii, Phys. Fluids 19 (1976) 1348. [8] R. Decoste and B.H. Ripin, Phys. Rev. Lett. 40 (1978) 34. [9] F.S. Felber and R. Decoste, Phys. Fluids 21(1978) 520. [10] A. Gurevich, D. Anderson and H. Wilhelmsson, Phys. Rev. Lett. 42 (1979) 769. [11] CE. Max, in: Laser-Plasma Interaction (North-Holland, Amsterdam, New York, Oxford. 1980) p. 301. [12] F. Begay and D.W. Forslund, Phys. Fluids 25 (1982) 1675. [13]P. Wagli and T.P. Donaldson, Phys. Rev. Lett. 40 (1978) 875. [14] P.M. Banks and T.E. Holzer, J. Geophys. Res. 73(1968) 6846. [15] N. Singh and R.W. Schunk, J. Geophys. Res. 87 (1982) 9154. [16] ED. Korop, Zh. Tekh. Fiz. 46 (1976) 2187 [Soy. Phys. Tech. Phys. 21(1976)1284]. [17] RB. Baksht, NA. Ratakin and BA. Kablamaev, Zh. Tekh. Fiz. 50 (1980) 487 [Soy. Phys. Tech. Phys. 25(1980) 294]. [18] B. Eiselt, Z. Phys. 132 (1952) 54. [19] IF. Kvarkhtsava, A.A. Plyutto, A.A. Chernov and V.V. Bondarenko, Soy. Phys. — JETP 3 (1956) 40. [20] R. Tanberg, Phys. Rev. 35(1930)1080. [21] A.A. Plyutto, Zh. Eksp. Teor. Fiz. 39 (1960) 1589 [Soy. Phys. — JETP 12 (1961) 1106]. [22] H.W. Hendel and T.T. Reboul, Phys. Fluids 5 (1962) 360. [23]P.F. Little, Journal of Nuclear Energy, Part C, 4 (1962) 15. [241A.A. Plyutto, V.N. Ryzhkov and AT. Kapin, Zh. Eksp. Teor. Fiz. 47 (1964) 494 [Soy. Phys. —JETP 20 (1965) 328]. [25] M.A. Tyulina, Zh. Tekh. Fiz. 35 (1965) 511 [Soy. Phys. Tech. Phys. 10 (1965) 396]. [26] P. Korn, T.C. Marshall and S.P. Schlesinger, Phys. Fluids 13 (1970) 517. [27] W. Demtröder and ‘A’. Jantz. Plasma Phys. 12 (1970) 691. [28] G.J. Tallents, Plasma Phys. 22 (1980) 709. [29] P. Pitsch, K. Rohr and H. Weber, J. Phys. D: AppI. Phys. 14 (1981) L51. [30] V.G. Eselevich and V.G. Fainshtain, Doki. Akad. Nauk SSSR 224 (1979) 1111 [Soy.Phys. Dokl. 24 (1979) 1141; Zh. Eksp. Teor. Fiz. 79 (1980) 870; [Soy. Phys. — JETP 52 (1980) 441]; Fiz. Plasmy 7 (1981) 503; [Soy.J. Plasma Phys. 7 (1981) 2711. [31] C. Chan, N. Hershkowitz, A. Ferreira, T. Intrator, B. Nelson and K. Lonngren, Phys. Fluids 27 (1984) 266. [32] A. Akhiezer, Zh. Eksp. Teor. Fiz. 47 (1964) 952 [Soy. Phys. — JETP 20 (1964) 637]. [33] A.V. Gurevich, LV. Pariiskaya and L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 49 (1965) 647 [Soy. Phys. — JETP 22 (1966) 449]. [34] D. Montgomery, Phys. Rev. Lett. 19 (1967) 1465. [35] A.V. Gurevich and L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 56 (1969) 1778 [Soy. Phys — JETP 29 (1969) 954]; Zh. Eksp. Teor. Fiz. 60(1971) 2155 [Soy. Phys. —JETP 33(1971)1159]. [36] J.E. Allen and J.G. Andrews, J. Phasma Phys. 4 (1970) 187.

394

Cli. Sack and H. Schamel, Plasma expansion into vacuum — A liydrodynamic approach

[37] M. Widner, J. Alexeff and W.D. iones, Phys. Fluids 14 (1971) 795. [38] R.J. Mason, Phys. Fluids 14 (1971) 1943. [39] PD. Prewett and J.E. Allen, I. Plasma Phys. 10 (1973) 451. [40] G.S. Voronov and L.E. Chernyshev, Zh. Tekh. Fiz. 43(1973)1484 [Soy.Phys. Tech. Phys. 18 (1974) 940]. [41]E.E. Lovetskii, AN. Polyanichev and VS. Fetisoy, Fiz. Plasmy 1 (1975) 773 [Soy.J. Plasma Phys. 1 (1975) 422]. [42]J.E. Crow, P.L. Auer and i.E. Allen, J. Plasma Phys. 14 (1975) 65. [431 A.V. Gurevich and L.P. Pitaeyskii, Prog. Aerospace Sci. 16 (1975) 227. [441 H. Shen and K.E. Lonngren, J. Engineering Math. 10 (1976) 135. [45] B. Bezzerides, D.W. Forslund and EL. Lindman, Phys. Fluids 21(1978) 2179. [46]J. Denavit, Phys. Fluids 22 (1979) 1384. [47] P. Mora and R. Pellat, Phys. Fluids 22 (1979) 2300. [48] K.E. Lonngren and N. Hershkowitz, IEEE Transactions on Plasma Science PS-7, No. 2 (1979) 107. [49]A.V. Gureyich and A.P. Meshcherkin, Zh. Eksp. Teor. Fiz. 53(1981)1810 [Soy. Phys. —JETP 53(1981) 937]. [50] MA. True, JR. Albritton and E.A. Williams, Phys. Fluids 24 (1981) 1885. [51] A.V. Gureyich and A.P. Meshcherkin, Fiz. Plasmy 8(1982)502 [Soy.J. Plasma Phys. 8(1982) 283]. [52] SI. Anisimoy, M.F. Ivanov, Yu.V. Medvedev and V.F. Shvets, Fiz. Plasmy 8 (1982) 1045 [Soy.I. Plasma Phys. 8 (1982) 591]. [53] L.D. Landau and EM. Lifschitz, Lehrbuch der theoretischen Physik Vol. VI (Akademie-Verlag, Berlin, 1966). [54] G.J. Barenblatt, Similarity, Self.Similarity, and Intermediate Asymptotics (Consultants Bureau, New York and London, 1979). [55] R.Z. Sagdeev. Reviews of Plasma Physics Vol. 4 (Consultants Bureau, New York, 1966) p. 23; J.E. Allen and L.M. Wickens. J. Physique, Colloque C7 (1979) 547. [56]Ch. Sack and H. Schamel, J. Comput. Phys. 53(1984) 395. [57]Ch. Sack and H. Schamel, Plasma Phys. Controlled Fusion 27 (1985) 717. [58] Ch. Sack and H. Schamel, Phys. Lctt. AllO (1985) 206. [59] R.L. Morse, Phys. Fluids 8 (1965) 308; D.W. Forslund and CR. Shonk, Phys. Rev. Lett. 25 (1970) 1699; H. Schamel, Plasma Phys. 14 (1972) 905. [60] W,F. Ames, Nonlinear Partial Differential Equations in Engineering, Mathematics in Science and Engineering Vol. 18 (Academic Press, New York, 1965). [61] K.E. Lonngren, in: Recent Advances in Plasma Physics, Proc. Workshop held at Physical Research Laboratory, Ahmedabad, ed. B. Buti (1977) p. 39. [62]R.C. Davidson, Methods in Nonlinear Plasma Theory (Academic Press, New York, 1972). [63] S.F. Johnson and K.E. Lonngren, Physica Scripta 25(1982)583. [64] H.L. Lewis, Preprint LA-UR-82-3100, Los Alamos National Laboratory, Los Alamos, New Mexico, USA, 1982. [65] RD. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, New York, 1967). [66] D. Potter, Computational Physics (Wiley, New York, 1973). [67] R.L. Dewar and E.J. Valeo, Proc. 6th Conf. on Num. Sim. of Plasmas, 1973. [68] D.W. Forsiund, J.M. Kindel, K. Lee and EL. Lindeman, Proc. 7th Conf. on Num. Sim. of Plasmas, Courant Institute NYU, 1975. [69] J. Trulsen, Tromsö Report, Tromsö, Norway, 1980. [70] A. Clebsch, J. Reine und Angew. Math. 56 (1859) 1. [71] K. Elsässer and H. Schamel, J. Plasma Phys. 15(1976) 409. [72] Ya. B. Zel’dovich and Yu.P. Raizer, Physics of Shock Waves and High Temperature Hydrodynamic Phenomena (Academic, New York, 1967); R. Courant and K.O. Friedrichs, Supersonic Flow and Shock Waves (Springer-Verlag, New York, Heidelberg, Berlin, 1976). [73] G. Kalman, Ann. of Physics 10 (1960) 1. [74] D.W. Forslund, i.M. Kindel, K. Lee and B.B. Godfrey, Phys. Fluids 22 (1979) 462. [75] N.S. Buchelnikova and E.P. Matochkin, Institute of Nuclear Physics, Report 83-88, Novosibirsk, USSR, 1983. [76] KY. Kim, Phys. Lett. A97 (1983) 45. [77] H. Schamel, Z. Naturforsch. 38a (1983) 1170. [78] P.K. Sakanaka, Phys. Fluids 15 (1972) 1323. [79] SI. Braginskii, Reviews of Plasma Physics Vol. 1 (Consultants Bureau, New York, 1965) p. 205. [80] A.J. Lichtenberg and MA. Lieberman, Regular and Stochastic Motion, Applied Mathematical Sciences Vol. 38 (Springer-Verlag, New York, Heidelberg, Berlin, 1983). [81] Ch. Sack, H. Schamel and R. Schmalz, Phys. Fluids 29 (1986) 1337. [82] E. Kamke, Differentialgleichungen, Losungsmethoden und Ldsungen I (Teubner.Verlag, Stuttgart, 1977). [83] B. Jüttner, XVII Intern. Conf. on Phenomena in Ionized Gases, Invited papers, Budapest (1985) p. 262. [84] V.A. lvanov, B. Jüttner and H. Pursch, Preprint 84-6, Akademie der Wissenschaften der DDR, Zentralinstitut für Elektronenphysik (1984). [85] GA. Mesyats, Proc. 10th Intern. Symp. Discharges and Electrical Insulation in Vacuum, Columbia (1982) p. 37. [86] DL. Book, NRL Plasma Formulary, Laboratory of Computational Physics, Naval Research Laboratory, Washington, D.C. 20375, 1980. [87] Yu.I. Chutov. XVI Intern. Conf. on Phenomena in Ionized Gases (ICPIG), Invited Papers, Düsseldorf, FRG, 1983, p303.

Ch. Sack and H. Schamel, Plasma expansion into vacuum — A hydrodynainic approach [88] L.P. Gorbachev, Magnetohydrodyn. 20 (1985) 398. [89]A.W. Ehler, J. AppI. Phys. 46 (1975) 2464. [901P.M. Campbell, R.R. Johnson, F.J. Mayer, LV. Powers and D.C. Slater, Phys. Rev. Let!. 39 (1977) 274. [91]1.5. Pearlman, J.J. Thomson and CE. Max, Phys. Rev. Lett. 38 (1977) 1397. [92]L.M. Wickens and J.E. Allen, Phys. Fluids 24 (1981) 1894. [93] G. Jordan-Engeln and F. Reutter, Numerische Mathematik für Ingenieure, BI-Htb 104, Mannheim, Wien, Zurich, FRG, 1973.

395