Two-dimensional particle simulation parallel electrodes
of plasma expansion
bletween plane
Kartik PateI=) and V. K. Mago Ion Extraction Section, B-204A Modular Labs, Laser and Plasma TechnologyDivision, Bhabha Atomic Research Centre, Trombay,Bombay 400-085, India
(Received 16 March 1995; acceptedfor publication 23 May 1995) We simulate in two dimensionsthe expansionof a plasma betweenbiasedplane parallel electrodes using the particle-in-ce&method. Such a plasma is frequently createdin many experimentsby the interaction of a pulsed laser with atomic vapor or gas stream. We describe the motion of the electrons and ions and reproducethe experimentally observedbulk drift of the plasma toward the high potential electrode.,This is explained in terms of the retrogrademovement of the ion sheath bound&y on one side of the plasma accompaniedby ambipolar diffusion on the other. We estimate the reducedion density in the plasmaby observing oscillations in the space-charge-limitedcurrent. By calculating the plasma decay time constant and ambipolar diffusion coefficient, we note that computersimulation can generatedata which is difficult to measureexperimentally,and not possible to calculate analytically. 0 1995 American Institute of Physics. 1. INTRODUCTION
Plasma simulation is widely used as a tool for carrying out plasmaphysics experimentson the computer.When done with care it can lead to significant savings in time and money. Not only can it generatedata similar to that which can be obtainedfrom actual experiments,it can also generate complementarydata which would -beimpossible or very difficult to acquire otherwise. A system consisting of an expandingplasmaboundedby plane parallel electrodesexists in an important class of experiments whose primary aim iS to understandthe spectral properties of atoms and molecules.The plasma in these experiments is formed usually by ionization of metal vapor or gas stream by single or multiple pulsed’laserbeams.Signals from theseexperimentsare dependenton the wavelengthand intensity of the incident radiation and on the collective behavior of the plasma in responseto external electromagnetic fields. In order to analyze such experimentsit becomesimportant, therefore, to understand the behavior of the plasma which is created in them. Such a plasma has characteristics which are different from plasmascreatedthermally or in discharge tubes. They are created almost instantaneously,are marked by steep density gradients,and have short lifetimes, ranging from 10 to 100 pus.The iodS usually have a directed velocity, since they are created within a beam, and all the electrons initially have the same energy since they are formed by ionization from a neutral speciesby incident laser radiation. In spacethe chargedistribution usually mirrors the intensity profile of the incident laser beams. This kind of plasma has been studied analytically by Chen,’Murakami and Nishihara,” and Okano,3 as well as experimentally by Giammanco,4 Ogura, Arisawa, and Shibata,5 Yamada, Tetsuka, and Deguchi,6 Hyman and Williamson,7 and Yamadaet aL8 Giammancohas also modeled this system using the fluid approximation. Elaborate computer simulations of this system have beencarried out in “Electronic mail:
[email protected]
J. Appt. Phys. 78 (7), 1 October 1995
one dimension by Watanabeand Okanogand in two dimensions by Ogura, Kaburaki, and Shibata,” both employing a hybrid approach.The fluid approachin one dimension has been used in simulations by Wid.neret al.” and in two dimensionsby Vitello, Cerjan, and Eraun.12A one-dimtnsional simulation using the particle model has been reported by Procassini, Birdsail, and Morsel3 and Calder, Hulbert, and Laframboise’4 have reported a particle simulation using a one-and-halfdimension code. While one-dimensional analysis and simulation faithfully reproducethe essential aspectsof the expansion, this approachmisses features which are essentially two dimensional in nature. We report the particle simulation in two dimensiohs of an expanding plasma created between plane parallel electrodesand its subsequentmovement under the influence of an externally applied electric field. This article sheds light on the processof removal of ions and explains some peculiar aspectsof the motion of the plasma in a direction against the electric field that have been observed experimentally4-6but not reported in a particle simulation. The article is organizedas follows: Section II describes the method of simulation, Sec. III describes the computational geometry and the plasma characteristics,and Sec. IV outlines the results obtained; these are subsequently discussedin Sec. V. II. PARTICLE-IN-CELL
PLASMA SIMULATION
The simulation was done using the particle-in-cell (PIC) technique in two dimensions.‘5-‘7 In this technique the plasma is approximatedby substituting the large number of plasma chargeswith a comparatively smaller number of particles. Each particle carries more mass and hence a proportionally greatercharge,with the chaxge-to-mass ratio remaining the same.The plasmais constructedout of theseparticles which move, subject to electromagnetic fields, within a simulation volume called the computational box. The computational box may, in principle, be of any shape, but for convenienceit is usually kept square,rectangular,or circular. The computationalEioxis discretizedinto a spatial mesh (or grid), by sets of lines parallel ‘to the boundaries.These
0021-8979/95/78(7)/4371/11/$6.00
0 1995 American institute of Physics
4371
Downloaded 19 Dec 2003 to 203.197.42.2. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/japo/japcr.jsp
TABLE I. Plasma parameters. Grotid side (0 volts)
Potential side (60 volts)
Charge density Electron temperature Ion temperature Ion drift velocity Plasma size Ion mass
l.OXIOy cmm3
3000 K 0 K (cold ions) 1 X lo5 cm s-l in positive y direction 10x 10 mm’137 (barium)
are solved to get the electric and magnetic fields at the mesh points. Third, thesevalues are usedto calculate the forces on the particles using the Lore&z equation
F=e(E+vxB). 4omm, x
Fourth, the Newton equations of motion are solved to find the new particle velocities
FIG. 1. Computational geometry showing the electrodes and initial plasma position. The length of the sides are 4 cm (horizontal, marked X) and 10 cm (vertical, marked Y). The left-hand-side electrode is held constant at 60 V. The right-hand-side electrode and the boundaries are kept at ground potential. The separation between the electrodes is 3.5 cm. The smaller square shows the position of the plasma of size 1 cm by 1 cm at time t=O. A-A’ is the line along which the potential is monitoied during thBsimulation (see Fig. 6).
lines are usually uniformly spaced.The intersection of these lines definesthe mesh points (sometimescalled grid points) on which the electromagnetic-fieldvariables (chargedensity, electric field, and magnetic-field intensity, etc.) are defined. The lines themselves define the cell boundaries, through which the current density is computed.The electromagnetic fields are found by numerically solving the discretized Maxwell equationson the mesh.At the start of the simulation the particles are initialized by assigning them positions and velocities in the desired manner.The particle positions and velocities are updatedat fixed intervals of time by solving the Lorentz force equation and the Newton equation of motion. The particles can thus move freely throughout the computational box. Interparticle binary collisions are neglected.Each particle interacts with others via the fields which are defined on the mesh points. Since there is no direct interaction among the particles in this model, it is known as the particle-mesh representation,and the method of simulation is called the+PICmethod. The computational procedure consists of the following sequentialstages.First, the initial positions and velocities of the particles are used to calculate the charge and current densities-onthe mesh points, using a linear weighing scheme for this purpose.Thesevalues are used as sourceterms in the field equations.Second, the field equations(in MKS units),
V.E=plq,,
(1)
V-B=O,
0) (3)
4372
J. Appl. Phys., Vol. 78, No. 7, 1 October 1995
(5)
dv F -=dt
(6>
m’
and finally the new velocities are used to updatethe particle positions which are carried forward into the next time step. The iterative execution of the above stagesresults in a selfconsistent evolution of the plasma. In a fully electromagnetic code, which solves for the self-consistent magnetic field and particle currents, all the four field equations would be required. Plasmas created in the mannerdescribedin Sec. I, however, can be treatedby a simpler system of equationsby introducing some approximations. These are (5)
ib)
(c)
neglecting effects dependenton AX/&, neglecting particle currents, andkeeping the magnetic field constant.
These approximationseliminate Eqs. (2), (3), and (4) above, with Poisson’sequationremaining as the only field equation. The system is reduced to the well-known electrostatic plasma.By exploiting symmetry along one of the dimensions the geometry can be reducedto two or even one dimension. This generally leads to a significant decreasein the computer resourcesrequired to solve the problem. The program used -in the current simulation solves the Poissonequationusing a five-point difference formula which is iterated using the successiveapproximation technique accompaniedby Chebyshevacceleration.This schemewas preferred over faster fast Fourier transform (FFQ-based meth-
TABLE II. Program input parameters. Electrode separation Electrode potential difference Magnetic field Number of particles Plasma frequency Time step (=OZ/plasma frequency)
35 mm 60 V OG 1 X lo5 negative and positive particles 1.78X109 rad s -I 1.12x10-~0 s
K. Pate1 and V. K. Mago
Downloaded 19 Dec 2003 to 203.197.42.2. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/japo/japcr.jsp
ods because of its flexibility in handling complex boundary conditions and presence of multiple electrodes within the computational box. III. COMPUTATIONAL
GEOMETRY
Figure 1 shows the rectangular computational box, of dimension 4X 10 cm”, used in the simulation. The mesh size is 40X100, that is it is divided into cells that are 1 mm in width. Vertical lines near the boundaries represent electrodes on which the external potentials are applied. They end 1 cm
above the bottom boundary of the box. The left-hand-side electrode is held at 60 V and the opposite electrode is kept at 0 V The boundaries are also kept at 0 V. The plasma is created with a slight offset toward the bottom of the box, as shown in Fig. 1. It is seeded with 2X105 particles, which are divided equally between the negative and positive particles. Each particle represents a large number of actual plasma charges,the exact value being calculated by the program at run time from the input data deffning the initial plasma size and charge density. Initially the particles are assigned random positions within the
id~~b~~~~iii::i;:~:I~~:~~
L*
-I
:
:
.,::
:.
~. \
:,
:,;
‘:
? .I
:
0 ‘:,
: : : ..
‘. i :
1:
,._ I : : : ;: ,. ., .
:::‘ :!
!*
ii
; ;
(ii :
‘::
:
.:.
.
A\
1: I. .
*
::
..a
:: ::
I :
I;;
j
.::
:
.
“’ ii;
:
;
.::
;
: “.
:: ::
,’
:i. ::
;
>’ .-’
i
.i
:
i ;
:‘:; :
0.
::
; J _.
.:.
-~
1 ,
-....m.. .
.
-.’
. ..
:
,’
. ,’
. ,
_..
:: ::
: :
. ( , i
:: ;.‘!
.
.*_.,,,,_...._..... f' Ill*lrllIr.iiiiiiijiii:;iiii;;!ti:!
iiifiii;iiIi’ii;iritiit;ir;ltr!ii’~:,
-".
.q-;.r:,
(4
. f
. .
. .
. .
. . . . :::::.::::
::...,
.I’
t -,
r. ‘,.31
r'llIIIIlI.i~~,*t~~t:r;Iflii;;i'!!i~~'
04
(d)
FIG. 2. Frames from the movie showing the position of the plasma at different times. The darker shade represents ions. The dotted lines are equipotentials. The first equipotential line (60 V) is coincident with the left-hand-side electrode. Moving toward the right-hand-side electrod’e, subsequent lines are at 30, 15, 6, 3, 1.5, 0.6 V, etc. This result in lines being bunched nearer to ground potential. The frames are at intervals of 20 000 time steps, corresponding to the following times: t= (a) 0, (b) 2.24, (c) 4.48, (d) 6.72, (e) 8.96, (f) 11.2, (g) 13.4, (h) 15.7, (i) 17.9, (j) 20.2, (k) 22.4, (1) 24.7, (m) 27.0, and (n) 29.2 ,u.s.
J. Appl. Phys., Vol. 78, No. 7, 1 October 1995
K. Pate1 and V. K. Mago
4373
Downloaded 19 Dec 2003 to 203.197.42.2. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/japo/japcr.jsp
. ..
:.,..~:::::‘.::...i-.“" ~rirrl~,,~~,:::III'iiiii;c:
-w>$I 9. 72; I*+, h
