Nonlinear dynamics in expanding plasmas

Volume 110A, number 4

PHYSICS LETTERS

22 July 1985

N O N L I N E A R D Y N A M I C S IN E X P A N D I N G P L A S M A S Ch. S A C K and H. S C H A M E L Institut fftr Theoretische Physik, Ruhr-Universiti~t Bochum, D-4630 Bochum 1, Fed. Rep. Germany

Received 5 October 1984; revised manuscript received 25 January 1985; accepted for publication 15 May 1985

The expansion of a plasma occupying initially a half-space is investigated numerically and, by means of a novel description of the ion fluid, also analytically. A simple wave structure is found in the collisionless approximation. Stabilized by dissipation, the associated ion bunching gives rise to a fast ion component, similar to the ion blow-off in laser fusion. Three non-stationary regimes of this strongest nonlinear and inhomogeneous dynamical system are distinguished and discussed. For large t the ion front propagates with a speed proportional to (t - t 1)1/2 where t 1 is a reference time. A simple picture emerges, explaining the diverse experimental data.

Plasma expansion is believed to play an important role in areas such as laser-matter interaction, pellet ablation as means of refuelling, plasma scrape-off in divertor or limiter tokamaks, arcing on surfaces, polar or stellar winds in astrophysics, etc. Belonging to one of the oldest disciplines in plasma physics, extensive numerical, analytical and experimental investigations have been carried out. In this letter we shall show first of all that an expanding plasma develops dynamical structures that are yet unexplored and not yet described in the plasma literature. In an earlier paper [1 ] we have described our hydrodynamic lagrangian expansion code and have interpreted a singularity occurring in the ion flow as being due to a numerical instability. Here, we stress on the mathematical and physical side o f the expansion problem and offer a new interpretation of this singular plasma behaviour, part of which has been published elsewhere [2]. We assume that the plasma occupies initially a halfspace, with some smooth, narrow transition, and that the ions are cold and without drift, initially. Normalized by the plasma quantities in the unperturbed region, the equations which govern the evolution, are the cold ion equations which are supplemented by Poisson's equation and by a Boltzmann law for the electrons: atn + ax(nU ) = O,

~2~b = n e - n, 206

~t v + OaxO = - a x ~ ,

ne(~ ) = exp(~b).

Fig. 1 is taken from ref. [1] and shows the spatial dependence of the ion density at different time steps; fig. la presents a global view whereas in fig. lb the leading front is drawn on a larger scale. We notice the appearance of a peaked ion front and the abovementioned singularity. A collapse of the ion density takes place at x = 4.11 and t = 17.94, which is called henceforth the critical point. Fig. 2 represents a local view of the front showing the final stage of the collapse in more detail. The ion velocity, plotted in fig. 2a, is seen to steepen progressively, and its negative gradient diverges at the critical point. The electric field, shown in fig. 2b behaves similarly with a reversed sign of the gradient. Interesting details of the ion density n, respectively o f its inverse, are found in fig. 2c where V = 1/n, the specific ion volume is displayed. It is V-shaped with two individual prongs, and its minimum is found to decrease linearly, reaching zero at the critical point. The electric potential and the electron density, not shown here, remain well behaved and smooth. Our first aim will be to show that this collapse reflects an intrinsic nonlinear process of the system, with other words, it is not of a numerical origin as believed earlier [1 ]. We do this by establishing for the first time a single scalar wave equation representing the entire system (1), and by solving this equation perturbatively near the critical point. We introduce the lagrangian mass variable

(1) 0.375-9601/85/$ 03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

Volume 110A, number 4

PHYSICS LETTERS

22 July 1985

n(x,t) 0.15~

n(x,t)

1.0

t=16 /t=12 ~t=O 075

1, \ \

\ \

t=13

o.,t \ \ \ ./

0.50

0.25 Ct

0

-2OO

.

.

.

.

i

-100

.

•

-

o

/-~

°.I

-20

-15

-10

-5

10

~X

-X

Fig. 1. Ion density as a function of space and time; (a) global behaviour of the density for three different time steps; the part framed by the rectangular is drawn on a larger scale, - 2 0 ~ x ~ 10, for five equidistant time steps in (b).

X

= f n(x', t) dx' = ~7(x, t),

(2)

0

-Z8 v(x,t)

-2.6

~

40o

~

~04

~p

x-b

from which, b y inversion, we get x = x(r/, t). Considering r / a n d t as the n e w i n d e p e n d e n t variables, the derivatives in eq. (1) t r a n s f o r m as a x ~ v - l an, a t + ua x ~ Ot, where V07, t) = 1/n(x07, t), t). F r o m this follows

V(rl, t) = Onx07 , t),

-Q22

f

/

[

/

Otx(~, t) = u(~, t ) ,

atU(rh t) = - V -10r~(rl , t) .

(3)

C o m b i n i n g these three equations we o b t a i n

(~==-a2t V = - a n v-lanq~. Poisson's e q u a t i o n , o n the other hand, becomes

016 I

I

I

X~

On V - 1 a ~

= neg-

1,

(5)

which, in comparison with eq. (4), yields

-70

ne(q~) = (1 - I7) V - 1 .

~0

(4)

__T~___k/____y_

(6)

Fig. 2. (a) Space-time dependence of the ion velocity o, (b) the electric field E and (c) the specific volume V = 1In, just before wavebreaking [taking place at (x., t.) = (4.112, 17.94)]. The minimum of the V-shaped specific volume approaches zero linearly (dashed line in (c)). 207

Volume ll0A, number 4

PHYSICS LETTERS

Inserting the Boltzmann law into eq. (6), the potential in eq. (4) can be replaced by ¢ = ln[(1 - 6 v - l ]

,

(7)

which finally gives 17+ a,7{(1 - 17)-1an [(1 - 17)/V]} = 0 .

(8)

We thus arrive at a formulation of the ion dynamical system (1) which, to our knowledge, has not been described and used before. This new scalar wave equation representing completely the ion dynamics inclusively charge separation is extremely nonlinear in contrast to the linear oscillator equation derived in a similar manner for finite amplitude electron Langmuir oscillations [3,4]. Hence, only approximative solutions of eq. (8) can be expected. As shown in detail in a subsequent paper [5], eq. (8) is satisfied by the following power series expansion V(r/, r) = a r [ l + r/2a - brl + (ne./6Fa 2) (*12 -- ~ 2 r 2 ) - (bne./18Fa 2 (~ - Far)2(r~ + 2Far) + O(e4)] ,

(9)

O <~ rl <~ g2r -- e ,

where a, b, Fa and ne, are constants, n e , representing the finite electron density at the critical point; r is given by r = t. - t, being our smallness quantity. Eq. (9) is an analytical expression for one branch of the V-shaped specific volume (see fig. 2c) in the r/, r space. It is valid near the critical point 77 = 0, r = 0. In order to get r/= 0 we have shifted the origin and replaced in eq. (2) x by ~ = x . - x. One easily verifies that the characteristic features of the dependent quantities of fig. 2 near the critical point are reproduced, which are V -+ 0 (resp. n -+ oo) a ~ v = - V -1 i)nv= - v - l

~nx = + ~ ' / V ~ + l / r ~ + ~ ,

= - V - 1 3nx = 1 7 / V ~ 1/ar -~oo,

O£E = - V - 1 3 n E

for r -+ 0+. In these expressions dot stands now for at, and E has changed sign, according to the abovecited transformation in x. Transforming back to the eulerian space which is accomplished by ~2r

~(rl, r ) = v , r -

we first get 208

f

drl' V ( n ' , r ) ,

22 July 1985

r / - Far = ~'(1 + b~/2 + hr) + O(e3), where h = bFa - (2a) -1 is constant and ~" := (Y - o . r ) / ar being of order e; o. is the value of o at the critical point. Insertion into eq. (9) yields V(~, r) = at{1 - hr - b~" + ~ [(ne./6922 ) (~ + 2Far) - b(b~/2 + hr)] } + O(e4),

(10)

which, replacing Y and r by x and t, represents one prong of the specific volume in the x - t space in the neighbourhood of the critical point. The other prong is given by another set of constants. Eq. (10) describes the spatial and temporal behaviour of V in fig. 2c. To lowest order it follows n(x, t) = [a(t. - t ) ] - I + 0 ( 1 ) , o(x, t) = ( x . - x ) / ( t . - t) + O(e) , E(x, t) = E . + O(e),

(11)

where again E . is the finite electric field at the critical point. The leading term in n is O(e - 1 ) and reflects ion bunching. To lowest order, namely O(e-1), the set of equations (11) satisfies iltV + Vi~xV = 0 .

(12)

This is the well-known equation for simple waves which intrinsically break. In the last stage of the collapse, therefore, the system behaves like an ordinary fluid, the electric field being dominated by the convective term oaxv which diverges. The numerical and analytical investigation, therefore, leads us to the conclusion that the dissipationless expansion of a plasma described by fluid equations is subject to ion wave-breaking. A similar conclusion was drawn for the quasi-neutral case by Akhiezer [6] in the cold ion limit assuming a homogeneous background plasma, and by Gurevich and Pitaevsky [7,8] for a finite ion temperature plasma expanding into vacuum. Ignoring for a moment the fact that quasi-neutrality breaks down at the steepened front, and using instead of eq. (8) the quasi-neutral version which becomes I7 + a 2 v - 1 = o, we obtain, of course, these quasi-neutral results, too. The singularity in the ion flow can thus be understood completely in terms of a simple wave structure and occurs whether quasi-neutrality is involved or not [5].

Volume 110A, number 4

PHYSICS LETTERS

Physically the collapse phenomenon can easily be understood. The ambipolar electric field which is located at the plasma front is hump-like and accelerated the ions producing in the early stage a similarly shaped velocity profile. This profile steepens at the front due to ion nonlinearity. Later on when the ion fluid has reached supersonic velocities, the convective term in the ion momentum equation overrules the electric field leading to a simple wave structure like in a usual hydrodynamic compressional wave. Let us next approach this problem from a kinetic point of view starting with a homogeneous cold ion plasma. Applying a sufficiently strong perturbation of the density or of the velocity, kinetic simulations [9,10] show the occurrence of ion wave-breaking and an associated ion bunching, and, as the time progresses, the formation of kinetic structures in the ion phase space. Depending on the strength of the initial perturbations these structures are found to be of double layer type for a weaker case (specifically, of the type of a slow ion acoustic double layer, in short SIADL [ 11 ]), of x-type for moderate strength and of multistreaming type for high strength [9,10]. As long as the ion velocity is single valued the plots of n, v and V arising from such simulations are seen to be locally identical with those we obtained. This observation we can evaluate in two ways. Firstly, we anticipate that a confident Vlasov simulation of the expansion problem will not differ from our numerical analysis up to the breaking point; furthermore, it will provide us an information about the state beyond wavebreaking, which in the collisionless regime has to be described in phase space rather than in real space. Secondly, the plasma-vacuum system can be understood as the extreme case of a density modulated (homogeneous) plasma, the density perturbation being so strong as to give zero density at the minimum. With other words, the plasmavacuum system is the strongest nonlinearly perturbed system we can imagine. Ion bunching and ion wave collapse are, therefore, natural events in such a system. Beyond wavebreaking the macroscopic description has to be given up, and a new world of kinetic phenomena is expected to occur which is not explored at all in this context. It should be noticed that particle simulations are probably not suited for resolving these structures due to the high noise level in the low density region where the steepening occurs.

22 July 1985

Fortunately, the abovementioned simulations also indicate a thermahzation of the ion distribution due to microinstabflities being generated by the structure. Thus, if we are interested in the long-term behaviour of the expansion, being aware that there is an intermediate stage where the evolution might not be described correctly, we can try to model these anomalous collision effects within a fluid description by adding a viscous term vO2v in the ion momentum equation [5]. A further argument supporting this model is the fact that for finite but small ion temperature the scale length in the front region can be much smaller than the mean free path of the ions, and hence multiple small angle Coulomb scattering will bring ion viscosity into play. Clearly, any hydrodynamic code which is able to go through this steepening and bunching phase must be dissipative as, e.g., Lax codes. To get an idea about the late phase of the expansion we have conducted several long-term calculations with a steeper initial density profde going beyond that of ref. [5]. Fig. 3 shows the space-time behaviour of the ion density for p = 5. We recognize that ion bunching is stabilized giving rise to a very pronounced ion peak. The viscosity has thus prevented the collapse. Identifying this peak with the component of fast ions observed in experiments, we conclude that the existence of these ions is a manifestation of the nonlinear processes taking place in the early stage of the plasma expansion. Near t = 95 (dashed line) the peaked ion density has reached a maximum value, and shortly after that it decreases significantly. Responsible for this delicate pattern is a bunching and debunching effect which can be understood by mass conservation and by the characteristic change of the velocity profile in this time interval. For t < 95 the velocity maximum Omax is behind the density spike implying the narrowing of the grid distances A in the ion peak region, and by mass conservation which in the Lagrangian picture reads n A = const, this means an increase of the local density. At t ~ 95 the velocity maximum overtakes the spike and lies for t > 95 in front of it. Consequently, A enlarges thereafter, and the density diminishes. We interpret this shift of Oma x as a tendency of the plasma to establish an overall self-similarity characterized by a velocity profde increasing linearly in x. It is, therefore, not dissipation which causes this rapid decrease of n but debunching. The instant where 209

Volume ll0A, number 4

PHYSICS LETTERS

22 July 1985

4O 2

1,oc~

~

'~

Qs-

/t gs

~

a6.

NR

0.4. 0.2.

0

w ....

0

-zoo

~0

~60

~0 ,

X

Fig. 3. Space-time evolution of the ion density in the viscid case, v = 5. Fez t = 95 (dashed curve) ion bunching is maximum, whereas for t > 95 debunching takes place.

Vma x passes the density peak, marks the onset o f a new, the final, acceleration mechanism giving rise to a further increase o f the fast ion velocity, as we shall discuss later. Coming back to fig. 3, we remark that a leading ion pulse is always there, although its value has be-

I-200 JL

.....

0I

i'

come small later on. Behind the front the scale lengths have widened, and the plasma has approached the quasi-neutral self-similar regime. This is visualized in fig. 4 showing the potential plot corresponding to fig. 3, in which this region is represented b y the straight portion o f ~. All the lines for different t have a tom-

i'200

i'

.oo

.....

T ' ~X

I

0

Fig. 4. Electrostatic potential as a function ofx for five different time steps. The positions of the ion front xf are indicated by arrows. For late times, the bulk plasma behaves quasi-neutral (straight lines of the potential curves). 210

Volume 110A, number 4

PHYSICS LETTERS E

v=5 .~

22 J u l y 1985

t = 175

• xD

I I

~10-2

-10

x1~-2

I

3 \ ' x ~ ss

/

"\'x

-8

v /'/"

~ - ~ "<"

E RAX

o=

f

o= o

2

%s _ - . . . .

.6

0.5

Es s

n

...

'~-

600

700

x

I ne

ox

• v = 0.2

-4

oV= = v=

1 5

-2

*

°x::.

lo

800

Fig. 5. I o n d e n s i t y n, e l e c t r o n d e n s i t y n e ( d o t t e d line), electric field E , a n d i o n v e l o c i t y o for v = 5 a t t = 175 ; o n l y the l e a d i n g f r o n t is d r a w n . F o r t h e sake o f c o m p a r i s o n , the corres p o n d i n g self-similar curves, d e n o t e d b y " s s ' , are p l o t t e d , too. T h e d o u b l e a r r o w a n d of r e p r e s e n t t h e l o c a l D e b y e l e n g t h a n d t h e i o n v e l o c i t y , r e s p e c t i v e l y , a t the p o s i t i o n o f m a x i m u m electric field.

It.

a

[:

5

10

20

:

:

; : ::.

50

: :

100

,I

I

I

200 300

t

10. xx

v,

.~'o°.. 8-

I I I

6-

mon fix point at ¢ = - l , x = - 40 (which is our front at t = 0) in accordance with the self-similar theory: ~b= - [(x + 40)/t + 1]. During the evolution the curvature of the potential which is initially positive corresponding to the pure electron cloud, reverses its overall sign in the front region indicating an ion rich region. This is in accordance with the laboratory experiments of Chan et al. [12] and indicates the violation of quasineutrality at the front for all times. To affirm this statement we examine the front further. Fig. 5 shows the spatial behaviour of the various quantities in this region. The self-similar solution corresponding to n, o and E, is denoted by nss, Oss and Ess. Its deviation from the true solution is salient. The electric field is found to be spatially nonconstant, the velocity increase is stronger than linear, and the densities, being smaller than nss , show the known characteristic deviations. The local Debye length given by n e 1/2, and indicated by an arrow at the position of maximum electric field, exceeds the local scale length of n. Consequently, quasi-neutrality and self-similarity do not apply to the front region and cannot be used to describe the front behaviour, not even in an approximative sense• On the other hand, as mentioned, the ion peak has become small, being only slightly separated from the bulk plasma. This leads us to the conclusion that in the last stage the dynamics of the major part of the plasma is determined by self-similarity. Only with

to :

• v=0.2 o V:

1

~ v:

5

I ~9"' "

I.¢ F°

[6 = ° "° "

TI'~"

I I

2

b 0 5

10

I I

I'/

......

I '°

20

S0

100

t

200 300

Fig. 6. (a) M a x i m u m electric field E m a x as a f u n c t i o n o f t i m e on a d o u b l e - l o g a r i t h m i c scale for t h r e e d i f f e r e n t values o f v, X : u = 5; o : u = 1 ; • : v = 0.2. (b) I o n f r o n t v e l o c i t y of as a f u n c t i o n o f t c o r r e s p o n d i n g t o (a). T h r e e r e g i m e s c a n b e d i s t i n g u i s h e d , s e p a r a t e d b y t . a n d t 0.

respect to the fastest ions the front behaviour has to be taken into account. We terminate the letter by investigating these fastest ions. In fig. 6 we plot the temporal behaviour of the maximum electric field E m a x (fig" 6a), and of the velocity of of the ion front (fig. 6b) defined by the local position ofEma x (see also fig. 5) for three values of v. Apparently three stages can be distinguished. The first one is given by the initial acceleration, steepening and bunching process in which the electric field is relatively large corresponding to the yet steep gradient. At t ~ t . , the point where in the absence of dissipation the collapse would take place, the second regime starts being characterized by a constant uf of approximately 3Cs, and a nearly constant maximum 211

Volume 110A, number 4

PHYSICS LETTERS

electric field. It is this stage where the fast ion peak is pronounced. Its duration depends on how fast the plasma establishes the self-similar motion. The process o f debunching sets in at the time t o where Vmax passes the ion peak. It terminates the second regime called, for short, the plateau regime. The final stage is characterized by a rapid increase of the front velocity and a decay of the maximum electric field Ema x. For large t the following formulae are found to be valid at the front [13]:

nf(t)=A/(tof(t) =

tl) ,

Emax(t)=

[8A(t- tl)]l/2 +B,

[2A/(t-

tl)] 1/2, (13)

where A, B, and t 1 are constants, slightly depending on v. For v = 0, we get by comparison with the numerical data,A = 8.76 × 1 0 - 3 , B = +6.17, and t 1 = 148.1. Eqs. (13) are derived by Taylor expansion of the full set of equations (1), at the front, x = xf(t), where the electric field is maximum, and by solving it under the assumption that E m ax(t) = (2 exp [q~f(t)] ) 1 / 2 [14] which is well satisfied in our numerical solution. Hence, the fast ion velocity increases asymptotically according to a square root law in the shifted time rather than to a logarithmic law as reported earlier, e.g. refs. [14,15] The effect o f the viscosity in fig. 6 is to increase the acceleration leading to a higher front velocity. Finally, we arrive at the following simple interpretation o f laboratory measurements, the results of which at a first glance, seem to be contradictory with respect to the maximum speed of the ions. There are, first o f all, experiments in which fast ion velocities of about 5c s and more were obtained [16,17]. Obviously, these experiments could reach the third regime which was made possible by the use o f a sufficiently long device (e.g. a length of 2 m). On the other hand, there is also a number of experiments [ 1 8 - 2 1 ] in which a maximum fast ion velocity of only 3c s was obtained. In view o f our numerical results it is clear how to interpret this limitation of of: The latter experiments namely, could only operate within the plateau regime. The reason is most probably the short distance (e.g. 2 0 - 5 0 era) and, hence, the limited time interval available for the expansion process. In conclusion, particular emphasis was given to the role ion nonlinearity plays during the evolution of an expanding plasma. The most salient feature in the hydrodynamic coUisionless regime was the occurrence o f ion wavebreaking and o f an ion density col212

22 July 1985

lapse indicating the transition to a kinetic ion phase space structure. Hydrodynamically, dissipation was able to stop the ion bunching creating a fast ion peak. Being, therefore, a relic of the nonlinear steepening process taking place at earlier times, the peak was found to be pronounced in an intermediate time interval only, the so-called plateau regime, and propagated with a velocity of approximately 3c s. For later times the fast ions were accelerated further, however, their number has decreased considerably, and the front has lost its importance in favour of an overall selfsimilar behaviour of the expanding plasma. The authors are grateful to Professor Dr. K. El~sser for valuable discussions. This work was supported by the Minister for Wissenschaft und Forschung des Landes Nordrhein-Westfalen under Grant No. IV B 4 - FA 9298.

[1] Ch. Sack and H. Schamel, J. Comput. Phys. 53 (1984) 395. [2] Ch. Sack and H. Schamel, Proc. Intern. Conf. on Plasma physics (Lausanne, 1984) pp. 4-6. [3] G. Kalman, Ann. Phys. (NY) 10 (1960) 1. [4] R.C. Davidson, Methods in nonlinear plasma theory (Academic Press, New York, 1972). [5] Ch. Sack and H. Schamel, Collapse of a plasma expanding into vacuum, preprint (1984). [6] I.A. Akhiezer, Soy. Phys. JETP 20 (1964) 637. [7] A.V. Gurevich and L.P. Pitaevskii, Soy. Phys. JETP 29 (1969) 954. [ 8 ] A.V. Gurevich and L.P. Pitaevskii, Prog. Aerospace Sci. 16 (1975) 227. [9] D.W. Forslund, J.M. Kindel, K. Lee and B.B. Godfrey, Phys. Fluids 22 (1979) 462. [10] N.S. Buehelnikova and E.P. Mathochkin, Institute of Nuclear Physics, Report 83-88, Novosibirsk (1983). [11] H. Scharnel, Z. Naturforsch. 38a (1983) 1170. [12] Ch. Chan, N. Hershkowitz, A. Ferreira, T. Intrator, B. Nelson and K. Lonngren, Phys. Fluids 27 (1984) 266. [13] Ch. Sack, dissertation, Ruhr-Universit//t Bochum (1985). [14] J.E. Crow, P.L. Auer and J.E. Allen, J. Plasma Phys. 14 (1975) 65. [15] J. Denavit, Phys. Fluids 22 (1979) 1384. [16] V.G. Eselevieh and V.G. Fainshtein, Soy. Phys. JETP 52 (1980) 441. [17] V.G. Eselevich and V.G. Fainshtein, Soy. J. Plasma Phys. 7 (1982) 271. [18] P.F. Little, J. Nucl. Energy C4 (1962) 15. [19] H.W. Hendel and T.T. Reboul, Phys. Fluids 5 (1962) 360. [20] A.A. Plyutto, V.N. Ryzhkov and A.T. Kapin, Soy. Phys. JETP 20 (1965) 328. [21] M.A. Tyulina, Soy. Phys. Tech. Phys. 10 (1965) 396.