Principles of âParticle in cellâ simulations
University of Ljubljana Faculty for mathematic and physics
Seminar
Principles of ”Particle in cell” simulations
Author: Gregor Filipič Mentor: dr. Tomaž Gyergyek April 2008
Abstract: In this seminar I will introduce the basics for computer particle plasma simulation. I will expose some problems we encounter when we simulate a large number of particles and how we can proceed. In the end I will introduce an example of plasma simulation and quick theoretical background of a problem.
1
Contents 1 Introduction
3
2 Computer simulations 2.1 Particle in cell simulations (PIC) . 2.2 Iteration steps in simulation . . . . 2.3 Integration of equations of motion 2.4 Integration of field equations . . . 2.5 Particle and force weighting . . . .
. . . . .
4 4 5 6 7 9
3 Example 3.1 Theoretical background - plasma sheath . . . . . . . . . . . . . . . . . . 3.2 Theoretical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Computer simulation result . . . . . . . . . . . . . . . . . . . . . . . . .
10 10 11 13
4 Conclusion
15
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
List of Figures 1 2 3 4 5 6
Computational cycle. Source [2] . . . . . . . . . . . . . . . . . . . . . . . Leap-frog integration sketch. Source [2] . . . . . . . . . . . . . . . . . . Idea of first- and second-order weighting and virtual shape of a particle they produce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Model of the plasma collector and source sheath in potential . . . . . . Resulting potential and particle densities from simulation . . . . . . . . Resulting particle distribution regarding velocity along x-axis for ions and electrons. In upper case is the real physical phenomena, and the lower graphs result from numerical effect. . . . . . . . . . . . . . . . . . . . . .
2
6 6 9 12 13
14
1
Introduction
Plasma is a system of large number of charged particles. They interact with each other through Coulomb force which is long ranged. That is why plasma shows collective behavior which is very interesting and physicists give great effort to understand them and possibility to foresee plasma behavior at specific conditions. We all know basic principles behind charged particle interaction and motion, they are described by Maxwell’s equations. It is easy to calculate how one charged particle will move if we switch on electric field, but for great number of particles is impossible to calculate the exact trajectories or response to outer disturbance. We can help ourself with introduction of distribution function, and then use equations of motion like Boltzmann’s or Vlasov’s (Boltzmann’s with some simplifications) 1 - which are approximate equation of motion. For better understanding we try to implicate basic physical lows and principles into physical model. Than we try to simulate that model with computer program. After seeing the result of simulation we can judge how appropriate the model is. Since physicists are quite skillful when it comes to computer programming there is no problem to compile physics into simulation program. But a problem is computer capability. We can not process infinite data in real time. We could simulate a few particles and get good results quickly. But number of particles in laboratory plasma is significantly larger, of order up to 1024 . The direct approach isn’t good. It is necessary to use some tricks and numerical computation. The first particle models of electrostatic plasma resulted from pioneering work of Buneman [7] at Standford on the ”dissipation of currents in ionized media” and of Birdsall and Bridges [8] at Berkley on ”space-charge instabilities in electron diodes and plasma converters.” These models were one dimensional and did not use mesh for calculation of the fields. First mesh models were independently developed at Standford by Burger(1965) in one dimension and Hockney (1996) in two dimension. Those algorithms used nearest grid point charge assignment and field interpolation in order to economize on computer time. Higher-order interpolation schemes were first used by the Berkley group (Birdsall and Fuss [9]) in order to reduce the noise in simulation. [1]
1
For details look my previous seminar: Uvodni pojmi iz fizike plazme ter nekaj zgledov iz magnetohidrodinamike
3
2
Computer simulations
There are two basic principles to start plasma simulations from: kinetic and fluid description. Kinetic simulations considers plasma modeling involving particle interactions through electromagnetic field. This is done by solving numerically kinetic plasma equations (Vlasov or Fokker-Planck equations) or with solving Maxwells equations for particles, their interactions with each other and interactions with externally applied electromagnetic field. It’s result is distribution function. This approach has been particularly successful in dealing with basic physical problems in which particle distributions deviate significantly from local Maxwellian. To reduce complexities of kinetic description, fluid model numerically solves the magnetohydrodynamic (MHD ) equations, which are obtained by taking velocity moments of Boltzmann equation. It describes plasma through macroscopic quantities (density, velocity, and mean energy). MHD principle is thus useful for large-scaled problems directly related to the behavior of experimental devices. For ”Particle in cell” simulations (PIC) we use kinetic principle. Plasma simulations started in the early 60’. In first step, it was naturally to simulate one dimensional plasma with relatively small number of particles. Since then the development of new algorithms and accessibility of more powerful computers has allowed particle simulations to progress to more complex and realistic situations with more particles and dimension. Especially great improvement in simulation capability was inventing parallel cluster computing (by Gene Amdahl, IBM in 1967, Amdahl’s Law - theoretical calculation of maximum expected improvement to an overall system, when only one part of the system is improved). [1] [2]
2.1
Particle in cell simulations (PIC)
The difference between real and simulated plasma lies in charge representation, the fields and the space-time in which phenomena occurs. In real plasma are usually electrons and ions which are charged positive or negative. In simulations we use computer particles. They represent fixed number of real particles. Thus they are heavier and have different charge than the real particles but they have the same charge to mass ratio. Number of particles in Debye cube ND = nλ3D can be from ≈ 102 in alkali vapor plasmas up to ≈ 106 in confined fusion plasmas. When we reduce this number in our simulation we still get good result if we are interested in collective behavior of collisionless plasma in system that is much longer than Debye length. In one dimensional simulation is required number of particles reduced to ND = nλ1D ≈ (106 )1/3 = 102 . We can alter the simulation so that we use even smaller ND but keep the essential plasma behavior. Motion of computer particles is determined in two steps. From initial currents and charge densities we calculate electromagnetic field by using Maxwell’s equations. Then we use this fields in Newton equation to move particles for a small distance accordingly. Then we recalculate the fields from new particle positions and charge density. We are
4
repeating this two steps during all simulation and simulate particle movement in mean field surrounding them. Rather than to solve equations of motion and Maxwellian equation (Coulomb’s law, Poisson’s equation) in continuous space and time, we divide the physical volume into cells by lines which run parallel to the boundaries. The intersections of these lines define set of points called mesh points or grid points. From this discretization comes the name for Particle in cell simulation. In these points we calculate the charge fields and relative to this points we move our particles. Their space coordinates are continuous and they can occupy position anywhere within the mesh. Local charge density is simply number of particles in individual cell and divided by volume of the same cell. Benefits of a spatial grid against direct Coulomb’s low are that we avoid the singularities at origin of each particles. Mesh points are also used in standard numerical solutions of differential equation. Besides, we are rarely interested in details of close encounters of particles and we don’t get electrical field by summing contribution of each particle but rather calculate local charge density, potential and then field. The mesh sets lower limit to spatial resolution of particle-particle and particle-field interaction. It is smart to keep length of individual spatial cell smaller then Debye length. Reason is, that Debye length is the distance over which charge density separation effect can occur, and if spatial cell is larger this effects are not seen in calculated fields. In other direction, we don’t want to have spatial cell shorter than the radius of computer particle, due to charge density calculation. [2] [3]
2.2
Iteration steps in simulation
Simulation works in steps. First step is to read the initial conditions, particle positions and velocities. In second step some weighting has to be done because calculation of charge density in mesh grid point depends on distance of particle from that point. After obtaining charge densities we proceed with integration of field equations. Next we weight how the field in individual grid point effect each particle. The last step is integration of equation of motion and to accordingly change position and velocity of particles. Then we repeat the procedure until we obtain the wanted result. Basic scheme of computational cycle is in Figure 1. In every step time changes for ∆t. The size of time step is chosen considering physical phenomena we are interested in. If we have small mass particles we will want to have shorter time step, because plasma frequency is bigger than with heavier particles for which we can afford longer time step. If we are particularly interested in low frequency effect we will use longer time step, which will shorten simulation time. In every cycle numerous integrations occurs. For every particle one integration of equation of motion and one field equation integration for each grid point. That gives for a problem that calls for 10000 particles to be run for 1000 time steps on a 1000 mesh grid points 1011 integrations. Thus choice of numerical method to use is very important. It has to be fast. Another thing to consider is computer memory. Minimum information needed for integration is particle velocity and position. In our problem that gives 20000 words. If we want to use high-order method (e.g., Runge-Kutta) we have to save records from several previous steps and the memory needed would increase, it would also multiply number of operations taken for each particle. That is why we almost always use the least information (storage) and the fastest method we can. [1] [2] 5
Figure 1: Computational cycle. Source [2]
2.3
Integration of equations of motion
One of commonly used is leap-frog method. We simply replace two first order differential equations with two finite-difference: m
d~v = F~ dt d~x = ~v dt
⇒ ⇒
~vnew − ~vold = F~old ∆t ~xnew − ~xold = ~vnew ∆t
m
(1) (2)
Because we need ~vnew in (eq. 1) from (eq. 2) we calculate position step at later time then velocity step (as sen in 2). We take this simulation flow into consider with initial velocity and position at t = 0. We push ~v (t = 0) back to ~v (− ∆t 2 ) using the force ~ F calculated at t = 0. Second energies calculated from v (kinetic) and x (potential, or field) must be adjusted to appear at the same time.
Figure 2: Leap-frog integration sketch. Source [2] The force F~ consists from two parts, ~ + q(~v × B) ~ F~ = F~electric + F~magnetic = q E 6
Magnetic and electric fields must be calculated at the particle. Therefore we must ~ and B ~ from the grid to the particle. interpolate E Let us consider an example of one dimensional simulation. Let there be displacement in x direction and have vx and vy and let it be uniform magnetic field B along z axis. Electric field E is directed in x direction. Magnetic force only rotate velocity while electric force changes amplitude of velocity in x direction. Velocity calculation by leap-frog method would look like this: • Half acceleration q ∆t ∆t vx (t′ ) = vx t − + Ex (t) 2 m 2 ∆t ′ vy (t ) = vy t − 2 • Rotation
vx (t′′ ) vy (t′′ )
=
cos ωc ∆t sin ωc ∆t − sin ωc ∆t cos ωc ∆t
• Half acceleration
∆t q ∆t vx t + = vx (t′′ ) + Ex (t) 2 m 2 ∆t = vy (t′′ ) vy t + 2 where t − ∆t/2 < t′ < t′′ < t + ∆t/2 and cyclotron frequency is ωc =
2.4
g m B.
Integration of field equations
From initial charge and current density we obtain electric and magnetic field using ~ = −∂ B/∂t ~ Maxwell’s equation. At this point we assume static case: ∇ × E ≈ 0 ~ = −∇Φ: following E ~ = −∇Φ ~ = ρ E ∇·E (3) ǫo Put this two together we obtain Poisson’s equation: ∇2 Φ = −
ρ ǫo
(4)
There are two approaches. First is to solve difference equations of (eq. 3) and (eq. 4): φj−1 − φj+1 2∆x ρj φj−1 − 2φj + φj+1 =− 2 (∆x) ǫo Ej =
7
(5) (6)
The later we can write in the matrix form as Aφ = −
(∆x)2 ρ ǫo
(7)
Using the boundary conditions and knowing all ρj we get as many equations as unknown φj . Second approach comes in handy when periodic conditions can be applied. In that case it is reasonable to use discrete Fourier transformation series for all grid quantities. ρ(~x) and Φ(~x) have Fourier transforms ρ(~k) and Φ(~k), where ~k is wave vector in the Fourier transform kernel, exp(−i~k · ~x). Thus second order derivate in Poisson’s equation is substituted by multiplying with −~k 2 , Φ(k) =
ρ(k) ǫo k 2
(8)
We have finite Fourier series; ρ(x) is known at positions of net grid points. Let denote net grid point with Xj , and the grid function (field, potential or charge density) with F (Xj ). Periodic condition for grid function is F (Xj ) = F (Xj + L). Now we can express finite Fourier transform as F (k) = ∆x
NX G−1
F (Xj )e−ikXj ,
(9)
j=0
where NG is total number of net grid points. The inverse Fourier transform is F (Xj ) =
1 L
N G/2−1
X
F (k)eikXj ,
(10)
n=−N G/2
Using above mathematical tool for finite difference equations (eq. 4) and (eq. 5), we obtain E(k) = −iκΦ(k) ρ(k) Φ(k) = ǫo K 2 where
sink∆x κ=k k∆x
2
and
K =k
2
sink∆x/2 k∆x/2
(11) (12) 2
Reducing the ∆x step size, κ approaches k and K 2 approaches k 2 . So finer grid means more accurate result. Basically we have one Fourier transformation of ρ(Xj ) → ρ(k), then one division by ǫo K 2 and multiplication with −iκ and the last step to obtain electrical field is inverse transform E(k) → E(Xj ). The wave vector k has values from kmin = 2π/L, where L is overall length of our system, to kmax = π/∆x and their negatives → λmin = 2∆x and λmax = L. [2]
8
2.5
Particle and force weighting
Through all simulation process we are calculating densities and fields at grid points from known positions of particles relative to grid points, and forces at particle position from nearest grid points field. This calculations are called weighting. That basically means that we are doing interpolation. This can be of a different order, depends on accuracy required. It is smart to use the same weighting for density and force calculation, otherwise self-force could happen and particles would accelerate them selfs. In zero-order weighting for density calculation at grid point noted by index j we simply count the number of particles (Nj ) in the grid point surroundings with radius ∆x/2 an divide it by volume of surrounding. In one dimensional simulation that means ∆x counting particles on interval [Xj − ∆x 2 , Xj + 2 ] and dividing with ∆x. Common name for this weighting is nearest grid point or NGP. Electric field used in the force is that at Xj for all the particles in the j th cell. It is quite fast because only one grid look-up is done. At the moment particles move out of j th cell at boundary Xj ± ∆x/2 density in j th cell falls for 1/∆x and by the same amount rises in (j + 1)th cell. It seems like the particle would have finite-size rectangular shape with length ∆x (Figure 3(a)). Density jumps at boundaries when particle passes through produce jumps in electric field also and both manifest in noise. This noise is intolerable in many plasma problems.
(a) Zero-order weighting
(b) First-order weighting
Figure 3: Idea of first- and second-order weighting and virtual shape of a particle they produce First order weighting smooths the density and field fluctuations, which reduces the noise. The down side is that uses more operations for the same job. It access two grid points for each particle, twice a step. The idea is, that particle contributes to j th and (j + 1)th cell charge density respectively to its position relative to Xj : ∆x − (xi − Xj ) Xj+1 − xi qj = qi = qi (13) ∆x ∆x xi − Xj (14) qj+1 = qi ∆x where xi and qi is a position of ith particle and its charge. Net effect is produced triangular shaped particle which has with 2∆x (Figure 3(b)). With linear interpolation for assignment particles to their nearest points would produce the same result. This 9
weighting is called particle-in-cell modeling. It gives resulting plasma density and field with much less noise than zero order weighting. If we would like to suppress noise even more, we should consider higher order, like quadratic and cubic splines but the cost of more computation is usually to high. The field or force weighting operates in the same manner. The zero-order NGP force comes from field at the nearest grid point. The first-order force PIC comes from linear interpolation equivalent to charge weighting: Xj+1 − xi xi − Xj Ej + Ej+1 (15) E(xj ) = ∆x ∆x Although first-order weighting consumes more computer time per particle it allows coarser grid and fewer particles than NGP for the same result. [1]
3
Example
We will simulate plasma between two electrodes. One (left) will be at floating potential and at the other will plasma particles flow into our system. First one is called collector and the other source electrode which potential we will set as zero. The velocity distribution will be Half-Maxwellian and we will use reflux 2 . We will be mostly interested in potential dependent of x-coordinate. We will use one-dimensional program: xpdp1, from The Plasma Theory and Simulation Group at Berkley [4].
3.1
Theoretical background - plasma sheath
As said plasma consists of charged particles, electrons and ions. This particles are light and quickly response to changes in electromagnetic field. Thus plasma is highly conducting, and we can consider it as equipotential at plasma potential φp . When plasma face an object (diagnostic probe or wall) spontaneous voltage difference develops. It is called floating potential φf and is usually negative relative to plasma potential. That happens due to mass and temperature difference between electrons and ions because usually electrons are faster (with higher temperature) and get to obstacle first and charge it negative. Slow ions than raise a potential for some degree while equilibrium establishes. The potential drop φf occurs in a thin layer-sheath between plasma and solid. The sheath thickness is of order the Debye length s ǫo kT , λD = e2o no where no is particle density. Plasma itself is quasi neutral ne = ni but the sheath has net positive charge per unit volume since the plasma electrons are repelled by negative potential on obstacle. The areal negative charge density in the sheath is approximately equal to positive charge density. The sheath thus shields the plasma from potential on the solid surface. This 2
that means that all the particles returned to the right electrodes will be returned in the system with the same velocity amplitude but opposite sign
10
shielding effect is important in probe plasma research when we have to consider that size of probe is effectively enlarged by plasma sheath. [5]
3.2
Theoretical model
Electrons are much lighter than ions thus they will reach collector first and charge it negative and it will start repelling the rest of electrons. They return to the source electrode. Slow varying potential drop occurs at the source at distance only few Debye lengths and forms a source sheath. After ions reach the collector equilibrium can be established. For calculating how potential looks we must solve Poisson’s equation but we have to know boundary conditions. We said that plane at x = 0 is reference potential of zero and no surface charge is allowed thus electric field is zero. Collector at x = L floats to potential Φc . Because collector is not connected to external world net electric current must be zero and the electron flux to collector must be equal to ion flux (eq. 19). The value of potential at the neutral or inflection point, where ∇2 Φ = 0, between the source and collector sheaths is designated Φp . The electric field −∇Φ at Φp , which by definition is a constant, is chosen to be zero when the source and collector sheaths are many Debye lengths apart. Thus with the zero source field, the total net charge enclosed from 0 ≤ Φ ≤ Φp , is zero. The two conditions of zero field at Φ = 0 and Φ = Φp , applied to the first integral of Poisson’s equation (eq. 20) and the neutrality condition at Φp (eq. 18), provide two equations with three unknowns: Φp , Φc and the ratio of emitted densities. [6] We model plasma with Maxwellian distribution functions for electrons and ions: r 1 me 2 exp − (1/2me v − eo φ) fe = noe 2πkTe kTe r mi 1 2 fi = noi exp − (1/2mi v + eo φ) (16) 2πkTi kTi We introduce variables mi eφ ψ= me kTe r v noi 2kTe u= vo = α= vo me noe and use them to write dimensionless distribution functions: p 1 Fe = √ exp −u2 + ψ H( ψ − ψc − u) π r p τµ e −µτ u2 − τ ψ H(− ψ/µ − u) Fi = α π τ=
Te Ti
µ=
(17)
With Heaviside function in eq. 17 we determine that ions travel only towards the collector and electrons can be repelled from the collector and also travel towards the source. Integrating distribution functions we get particle densities and zero charge condition at plasma potential φp : p p ψp − ψc = αe−τ ψp Erf c −τ ψp (18) eψp 1 + Erf 11
For zero net flux at collector we need next higher velocity momentum and obtain simple equation: α 1 (19) − √ eψc = − √ 2 πµτ 2 π Dimensionless Poisson equation and zero-field condition give us third equation: p p 1 2 p 2 p ψp ψc ψc −ψp √ √ 1−e − −ψc e + ψp − ψc e + Erf −ψc − e ψp − ψc = Erf 2 π π r p τ α 1 − e−τ ψp − 2 − ψp + e−τ ψp Erf −τ ψp (20) =− 2τ π
To obtain values of ψp , ψc and α for given τ and µ we numerically solve system of equations (18-20). The resulting potential distribution shape is shown on next scheme:
Figure 4: Model of the plasma collector and source sheath in potential
Numerical solution for µ = 5 ∗ 10−4 gives us: τ 10 100 1000
α 1,25 1,81 2,38
ψp -0,09 -0,21 -0,25
ψc -0,62 -0,69 -0,74
Table 1: Numerical results
12
3.3
Computer simulation result
When we simulate the same problem in xpdp1 we get the following results:
(a) Potential
(b) Particle density
Figure 5: Resulting potential and particle densities from simulation We see that theoretical model and plasma simulation for bounded plasma agree. The last step in the research is to measure plasma potential in laboratory and if data suits above mentioned methods, we can be pretty sure we know what is going on in plasma in our system.
13
In lower figure 6 is particle distribution through our system depending on their velocities. In figure 6(a) and 6(b) are shown good simulation and in figure 6(c) and 6(d) are results of inappropriate chosen parameters - if we choose too large time step or not enough particle cells for given system length (or possibly both).
(a) Electron velocity
(b) Ion velocity
(c) Electron velocity with bad param- (d) Ion velocity with bad parameters eters
Figure 6: Resulting particle distribution regarding velocity along x-axis for ions and electrons. In upper case is the real physical phenomena, and the lower graphs result from numerical effect.
14
4
Conclusion
I introduced basic one dimensional computer simulation in plasma physics. We saw, how particle flux equilibrium is achieved, although due to lack of collisions, we don’t have any thermodynamical relaxation (besides artificial re-Maxwellization through reflux). At this time more realistic simulations are made and constantly upgraded. Simulation nowadays are not limited to one dimension, interaction between particles is not only through electromagnetic field but also by direct collisions. Various cross sections for collisions are determined and integrated into simulations. It is also taken into account collisions of charged particles with neutrals, and ionization within plasma. Next step is to consider secondary emission after particle collide into electrode and knocks out an electron. To shorten simulation time and improve accuracy variable grids are introduced - with fewer points at constant potential and more points at stirred up potential. To speed up simulation and capability to simulate more particles, programs are adjusted to run simultaneous on several computer processors (clusters). New numerical methods can also improve the simulation performance.
15
References [1] R.W.Hockney, J.W.Eastwood, Computer simulation using particles (IOP Publishing Ltd, 1999) [2] C.K.Birdsall, A.B.Langdon, Plasma physics via computer simulation (McGrrawHill, Inc., 1985) [3] http://www.geocities.com/letapk/pic.html , april 2008 [4] http://ptsg.eecs.berkeley.edu/ , april 2008 [5] P.C.Stangeby, The Plasma Boundary of Magnetic Fusion Devices (Taylor & Francis, 2000) [6] L.A.Schwager, C.K.Birdsall, Collector and source sheats of finite ion temperature plasma (Physics of Fluids B: Plasma Physics, 1990) [7] O.Buneman, Dissipation of Currents in ionized media (Phys. Rev. 115:503, 1959) [8] C.K.Birdsall, W.B.Bridges, Space-charge instabilities in electron diodes and plasma converters (J. Applied Physics, vol. 32, no. 12, pp. 2611-2618, Dec. 1961) [9] C.K.Birdsall, D.Fuss, Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation (J. Comput. Phys. 3, 494–511, 1969)
16
Seminar
Principles of ”Particle in cell” simulations
Author: Gregor Filipič Mentor: dr. Tomaž Gyergyek April 2008
Abstract: In this seminar I will introduce the basics for computer particle plasma simulation. I will expose some problems we encounter when we simulate a large number of particles and how we can proceed. In the end I will introduce an example of plasma simulation and quick theoretical background of a problem.
1
Contents 1 Introduction
3
2 Computer simulations 2.1 Particle in cell simulations (PIC) . 2.2 Iteration steps in simulation . . . . 2.3 Integration of equations of motion 2.4 Integration of field equations . . . 2.5 Particle and force weighting . . . .
. . . . .
4 4 5 6 7 9
3 Example 3.1 Theoretical background - plasma sheath . . . . . . . . . . . . . . . . . . 3.2 Theoretical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Computer simulation result . . . . . . . . . . . . . . . . . . . . . . . . .
10 10 11 13
4 Conclusion
15
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
List of Figures 1 2 3 4 5 6
Computational cycle. Source [2] . . . . . . . . . . . . . . . . . . . . . . . Leap-frog integration sketch. Source [2] . . . . . . . . . . . . . . . . . . Idea of first- and second-order weighting and virtual shape of a particle they produce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Model of the plasma collector and source sheath in potential . . . . . . Resulting potential and particle densities from simulation . . . . . . . . Resulting particle distribution regarding velocity along x-axis for ions and electrons. In upper case is the real physical phenomena, and the lower graphs result from numerical effect. . . . . . . . . . . . . . . . . . . . . .
2
6 6 9 12 13
14
1
Introduction
Plasma is a system of large number of charged particles. They interact with each other through Coulomb force which is long ranged. That is why plasma shows collective behavior which is very interesting and physicists give great effort to understand them and possibility to foresee plasma behavior at specific conditions. We all know basic principles behind charged particle interaction and motion, they are described by Maxwell’s equations. It is easy to calculate how one charged particle will move if we switch on electric field, but for great number of particles is impossible to calculate the exact trajectories or response to outer disturbance. We can help ourself with introduction of distribution function, and then use equations of motion like Boltzmann’s or Vlasov’s (Boltzmann’s with some simplifications) 1 - which are approximate equation of motion. For better understanding we try to implicate basic physical lows and principles into physical model. Than we try to simulate that model with computer program. After seeing the result of simulation we can judge how appropriate the model is. Since physicists are quite skillful when it comes to computer programming there is no problem to compile physics into simulation program. But a problem is computer capability. We can not process infinite data in real time. We could simulate a few particles and get good results quickly. But number of particles in laboratory plasma is significantly larger, of order up to 1024 . The direct approach isn’t good. It is necessary to use some tricks and numerical computation. The first particle models of electrostatic plasma resulted from pioneering work of Buneman [7] at Standford on the ”dissipation of currents in ionized media” and of Birdsall and Bridges [8] at Berkley on ”space-charge instabilities in electron diodes and plasma converters.” These models were one dimensional and did not use mesh for calculation of the fields. First mesh models were independently developed at Standford by Burger(1965) in one dimension and Hockney (1996) in two dimension. Those algorithms used nearest grid point charge assignment and field interpolation in order to economize on computer time. Higher-order interpolation schemes were first used by the Berkley group (Birdsall and Fuss [9]) in order to reduce the noise in simulation. [1]
1
For details look my previous seminar: Uvodni pojmi iz fizike plazme ter nekaj zgledov iz magnetohidrodinamike
3
2
Computer simulations
There are two basic principles to start plasma simulations from: kinetic and fluid description. Kinetic simulations considers plasma modeling involving particle interactions through electromagnetic field. This is done by solving numerically kinetic plasma equations (Vlasov or Fokker-Planck equations) or with solving Maxwells equations for particles, their interactions with each other and interactions with externally applied electromagnetic field. It’s result is distribution function. This approach has been particularly successful in dealing with basic physical problems in which particle distributions deviate significantly from local Maxwellian. To reduce complexities of kinetic description, fluid model numerically solves the magnetohydrodynamic (MHD ) equations, which are obtained by taking velocity moments of Boltzmann equation. It describes plasma through macroscopic quantities (density, velocity, and mean energy). MHD principle is thus useful for large-scaled problems directly related to the behavior of experimental devices. For ”Particle in cell” simulations (PIC) we use kinetic principle. Plasma simulations started in the early 60’. In first step, it was naturally to simulate one dimensional plasma with relatively small number of particles. Since then the development of new algorithms and accessibility of more powerful computers has allowed particle simulations to progress to more complex and realistic situations with more particles and dimension. Especially great improvement in simulation capability was inventing parallel cluster computing (by Gene Amdahl, IBM in 1967, Amdahl’s Law - theoretical calculation of maximum expected improvement to an overall system, when only one part of the system is improved). [1] [2]
2.1
Particle in cell simulations (PIC)
The difference between real and simulated plasma lies in charge representation, the fields and the space-time in which phenomena occurs. In real plasma are usually electrons and ions which are charged positive or negative. In simulations we use computer particles. They represent fixed number of real particles. Thus they are heavier and have different charge than the real particles but they have the same charge to mass ratio. Number of particles in Debye cube ND = nλ3D can be from ≈ 102 in alkali vapor plasmas up to ≈ 106 in confined fusion plasmas. When we reduce this number in our simulation we still get good result if we are interested in collective behavior of collisionless plasma in system that is much longer than Debye length. In one dimensional simulation is required number of particles reduced to ND = nλ1D ≈ (106 )1/3 = 102 . We can alter the simulation so that we use even smaller ND but keep the essential plasma behavior. Motion of computer particles is determined in two steps. From initial currents and charge densities we calculate electromagnetic field by using Maxwell’s equations. Then we use this fields in Newton equation to move particles for a small distance accordingly. Then we recalculate the fields from new particle positions and charge density. We are
4
repeating this two steps during all simulation and simulate particle movement in mean field surrounding them. Rather than to solve equations of motion and Maxwellian equation (Coulomb’s law, Poisson’s equation) in continuous space and time, we divide the physical volume into cells by lines which run parallel to the boundaries. The intersections of these lines define set of points called mesh points or grid points. From this discretization comes the name for Particle in cell simulation. In these points we calculate the charge fields and relative to this points we move our particles. Their space coordinates are continuous and they can occupy position anywhere within the mesh. Local charge density is simply number of particles in individual cell and divided by volume of the same cell. Benefits of a spatial grid against direct Coulomb’s low are that we avoid the singularities at origin of each particles. Mesh points are also used in standard numerical solutions of differential equation. Besides, we are rarely interested in details of close encounters of particles and we don’t get electrical field by summing contribution of each particle but rather calculate local charge density, potential and then field. The mesh sets lower limit to spatial resolution of particle-particle and particle-field interaction. It is smart to keep length of individual spatial cell smaller then Debye length. Reason is, that Debye length is the distance over which charge density separation effect can occur, and if spatial cell is larger this effects are not seen in calculated fields. In other direction, we don’t want to have spatial cell shorter than the radius of computer particle, due to charge density calculation. [2] [3]
2.2
Iteration steps in simulation
Simulation works in steps. First step is to read the initial conditions, particle positions and velocities. In second step some weighting has to be done because calculation of charge density in mesh grid point depends on distance of particle from that point. After obtaining charge densities we proceed with integration of field equations. Next we weight how the field in individual grid point effect each particle. The last step is integration of equation of motion and to accordingly change position and velocity of particles. Then we repeat the procedure until we obtain the wanted result. Basic scheme of computational cycle is in Figure 1. In every step time changes for ∆t. The size of time step is chosen considering physical phenomena we are interested in. If we have small mass particles we will want to have shorter time step, because plasma frequency is bigger than with heavier particles for which we can afford longer time step. If we are particularly interested in low frequency effect we will use longer time step, which will shorten simulation time. In every cycle numerous integrations occurs. For every particle one integration of equation of motion and one field equation integration for each grid point. That gives for a problem that calls for 10000 particles to be run for 1000 time steps on a 1000 mesh grid points 1011 integrations. Thus choice of numerical method to use is very important. It has to be fast. Another thing to consider is computer memory. Minimum information needed for integration is particle velocity and position. In our problem that gives 20000 words. If we want to use high-order method (e.g., Runge-Kutta) we have to save records from several previous steps and the memory needed would increase, it would also multiply number of operations taken for each particle. That is why we almost always use the least information (storage) and the fastest method we can. [1] [2] 5
Figure 1: Computational cycle. Source [2]
2.3
Integration of equations of motion
One of commonly used is leap-frog method. We simply replace two first order differential equations with two finite-difference: m
d~v = F~ dt d~x = ~v dt
⇒ ⇒
~vnew − ~vold = F~old ∆t ~xnew − ~xold = ~vnew ∆t
m
(1) (2)
Because we need ~vnew in (eq. 1) from (eq. 2) we calculate position step at later time then velocity step (as sen in 2). We take this simulation flow into consider with initial velocity and position at t = 0. We push ~v (t = 0) back to ~v (− ∆t 2 ) using the force ~ F calculated at t = 0. Second energies calculated from v (kinetic) and x (potential, or field) must be adjusted to appear at the same time.
Figure 2: Leap-frog integration sketch. Source [2] The force F~ consists from two parts, ~ + q(~v × B) ~ F~ = F~electric + F~magnetic = q E 6
Magnetic and electric fields must be calculated at the particle. Therefore we must ~ and B ~ from the grid to the particle. interpolate E Let us consider an example of one dimensional simulation. Let there be displacement in x direction and have vx and vy and let it be uniform magnetic field B along z axis. Electric field E is directed in x direction. Magnetic force only rotate velocity while electric force changes amplitude of velocity in x direction. Velocity calculation by leap-frog method would look like this: • Half acceleration q ∆t ∆t vx (t′ ) = vx t − + Ex (t) 2 m 2 ∆t ′ vy (t ) = vy t − 2 • Rotation
vx (t′′ ) vy (t′′ )
=
cos ωc ∆t sin ωc ∆t − sin ωc ∆t cos ωc ∆t
• Half acceleration
∆t q ∆t vx t + = vx (t′′ ) + Ex (t) 2 m 2 ∆t = vy (t′′ ) vy t + 2 where t − ∆t/2 < t′ < t′′ < t + ∆t/2 and cyclotron frequency is ωc =
2.4
g m B.
Integration of field equations
From initial charge and current density we obtain electric and magnetic field using ~ = −∂ B/∂t ~ Maxwell’s equation. At this point we assume static case: ∇ × E ≈ 0 ~ = −∇Φ: following E ~ = −∇Φ ~ = ρ E ∇·E (3) ǫo Put this two together we obtain Poisson’s equation: ∇2 Φ = −
ρ ǫo
(4)
There are two approaches. First is to solve difference equations of (eq. 3) and (eq. 4): φj−1 − φj+1 2∆x ρj φj−1 − 2φj + φj+1 =− 2 (∆x) ǫo Ej =
7
(5) (6)
The later we can write in the matrix form as Aφ = −
(∆x)2 ρ ǫo
(7)
Using the boundary conditions and knowing all ρj we get as many equations as unknown φj . Second approach comes in handy when periodic conditions can be applied. In that case it is reasonable to use discrete Fourier transformation series for all grid quantities. ρ(~x) and Φ(~x) have Fourier transforms ρ(~k) and Φ(~k), where ~k is wave vector in the Fourier transform kernel, exp(−i~k · ~x). Thus second order derivate in Poisson’s equation is substituted by multiplying with −~k 2 , Φ(k) =
ρ(k) ǫo k 2
(8)
We have finite Fourier series; ρ(x) is known at positions of net grid points. Let denote net grid point with Xj , and the grid function (field, potential or charge density) with F (Xj ). Periodic condition for grid function is F (Xj ) = F (Xj + L). Now we can express finite Fourier transform as F (k) = ∆x
NX G−1
F (Xj )e−ikXj ,
(9)
j=0
where NG is total number of net grid points. The inverse Fourier transform is F (Xj ) =
1 L
N G/2−1
X
F (k)eikXj ,
(10)
n=−N G/2
Using above mathematical tool for finite difference equations (eq. 4) and (eq. 5), we obtain E(k) = −iκΦ(k) ρ(k) Φ(k) = ǫo K 2 where
sink∆x κ=k k∆x
2
and
K =k
2
sink∆x/2 k∆x/2
(11) (12) 2
Reducing the ∆x step size, κ approaches k and K 2 approaches k 2 . So finer grid means more accurate result. Basically we have one Fourier transformation of ρ(Xj ) → ρ(k), then one division by ǫo K 2 and multiplication with −iκ and the last step to obtain electrical field is inverse transform E(k) → E(Xj ). The wave vector k has values from kmin = 2π/L, where L is overall length of our system, to kmax = π/∆x and their negatives → λmin = 2∆x and λmax = L. [2]
8
2.5
Particle and force weighting
Through all simulation process we are calculating densities and fields at grid points from known positions of particles relative to grid points, and forces at particle position from nearest grid points field. This calculations are called weighting. That basically means that we are doing interpolation. This can be of a different order, depends on accuracy required. It is smart to use the same weighting for density and force calculation, otherwise self-force could happen and particles would accelerate them selfs. In zero-order weighting for density calculation at grid point noted by index j we simply count the number of particles (Nj ) in the grid point surroundings with radius ∆x/2 an divide it by volume of surrounding. In one dimensional simulation that means ∆x counting particles on interval [Xj − ∆x 2 , Xj + 2 ] and dividing with ∆x. Common name for this weighting is nearest grid point or NGP. Electric field used in the force is that at Xj for all the particles in the j th cell. It is quite fast because only one grid look-up is done. At the moment particles move out of j th cell at boundary Xj ± ∆x/2 density in j th cell falls for 1/∆x and by the same amount rises in (j + 1)th cell. It seems like the particle would have finite-size rectangular shape with length ∆x (Figure 3(a)). Density jumps at boundaries when particle passes through produce jumps in electric field also and both manifest in noise. This noise is intolerable in many plasma problems.
(a) Zero-order weighting
(b) First-order weighting
Figure 3: Idea of first- and second-order weighting and virtual shape of a particle they produce First order weighting smooths the density and field fluctuations, which reduces the noise. The down side is that uses more operations for the same job. It access two grid points for each particle, twice a step. The idea is, that particle contributes to j th and (j + 1)th cell charge density respectively to its position relative to Xj : ∆x − (xi − Xj ) Xj+1 − xi qj = qi = qi (13) ∆x ∆x xi − Xj (14) qj+1 = qi ∆x where xi and qi is a position of ith particle and its charge. Net effect is produced triangular shaped particle which has with 2∆x (Figure 3(b)). With linear interpolation for assignment particles to their nearest points would produce the same result. This 9
weighting is called particle-in-cell modeling. It gives resulting plasma density and field with much less noise than zero order weighting. If we would like to suppress noise even more, we should consider higher order, like quadratic and cubic splines but the cost of more computation is usually to high. The field or force weighting operates in the same manner. The zero-order NGP force comes from field at the nearest grid point. The first-order force PIC comes from linear interpolation equivalent to charge weighting: Xj+1 − xi xi − Xj Ej + Ej+1 (15) E(xj ) = ∆x ∆x Although first-order weighting consumes more computer time per particle it allows coarser grid and fewer particles than NGP for the same result. [1]
3
Example
We will simulate plasma between two electrodes. One (left) will be at floating potential and at the other will plasma particles flow into our system. First one is called collector and the other source electrode which potential we will set as zero. The velocity distribution will be Half-Maxwellian and we will use reflux 2 . We will be mostly interested in potential dependent of x-coordinate. We will use one-dimensional program: xpdp1, from The Plasma Theory and Simulation Group at Berkley [4].
3.1
Theoretical background - plasma sheath
As said plasma consists of charged particles, electrons and ions. This particles are light and quickly response to changes in electromagnetic field. Thus plasma is highly conducting, and we can consider it as equipotential at plasma potential φp . When plasma face an object (diagnostic probe or wall) spontaneous voltage difference develops. It is called floating potential φf and is usually negative relative to plasma potential. That happens due to mass and temperature difference between electrons and ions because usually electrons are faster (with higher temperature) and get to obstacle first and charge it negative. Slow ions than raise a potential for some degree while equilibrium establishes. The potential drop φf occurs in a thin layer-sheath between plasma and solid. The sheath thickness is of order the Debye length s ǫo kT , λD = e2o no where no is particle density. Plasma itself is quasi neutral ne = ni but the sheath has net positive charge per unit volume since the plasma electrons are repelled by negative potential on obstacle. The areal negative charge density in the sheath is approximately equal to positive charge density. The sheath thus shields the plasma from potential on the solid surface. This 2
that means that all the particles returned to the right electrodes will be returned in the system with the same velocity amplitude but opposite sign
10
shielding effect is important in probe plasma research when we have to consider that size of probe is effectively enlarged by plasma sheath. [5]
3.2
Theoretical model
Electrons are much lighter than ions thus they will reach collector first and charge it negative and it will start repelling the rest of electrons. They return to the source electrode. Slow varying potential drop occurs at the source at distance only few Debye lengths and forms a source sheath. After ions reach the collector equilibrium can be established. For calculating how potential looks we must solve Poisson’s equation but we have to know boundary conditions. We said that plane at x = 0 is reference potential of zero and no surface charge is allowed thus electric field is zero. Collector at x = L floats to potential Φc . Because collector is not connected to external world net electric current must be zero and the electron flux to collector must be equal to ion flux (eq. 19). The value of potential at the neutral or inflection point, where ∇2 Φ = 0, between the source and collector sheaths is designated Φp . The electric field −∇Φ at Φp , which by definition is a constant, is chosen to be zero when the source and collector sheaths are many Debye lengths apart. Thus with the zero source field, the total net charge enclosed from 0 ≤ Φ ≤ Φp , is zero. The two conditions of zero field at Φ = 0 and Φ = Φp , applied to the first integral of Poisson’s equation (eq. 20) and the neutrality condition at Φp (eq. 18), provide two equations with three unknowns: Φp , Φc and the ratio of emitted densities. [6] We model plasma with Maxwellian distribution functions for electrons and ions: r 1 me 2 exp − (1/2me v − eo φ) fe = noe 2πkTe kTe r mi 1 2 fi = noi exp − (1/2mi v + eo φ) (16) 2πkTi kTi We introduce variables mi eφ ψ= me kTe r v noi 2kTe u= vo = α= vo me noe and use them to write dimensionless distribution functions: p 1 Fe = √ exp −u2 + ψ H( ψ − ψc − u) π r p τµ e −µτ u2 − τ ψ H(− ψ/µ − u) Fi = α π τ=
Te Ti
µ=
(17)
With Heaviside function in eq. 17 we determine that ions travel only towards the collector and electrons can be repelled from the collector and also travel towards the source. Integrating distribution functions we get particle densities and zero charge condition at plasma potential φp : p p ψp − ψc = αe−τ ψp Erf c −τ ψp (18) eψp 1 + Erf 11
For zero net flux at collector we need next higher velocity momentum and obtain simple equation: α 1 (19) − √ eψc = − √ 2 πµτ 2 π Dimensionless Poisson equation and zero-field condition give us third equation: p p 1 2 p 2 p ψp ψc ψc −ψp √ √ 1−e − −ψc e + ψp − ψc e + Erf −ψc − e ψp − ψc = Erf 2 π π r p τ α 1 − e−τ ψp − 2 − ψp + e−τ ψp Erf −τ ψp (20) =− 2τ π
To obtain values of ψp , ψc and α for given τ and µ we numerically solve system of equations (18-20). The resulting potential distribution shape is shown on next scheme:
Figure 4: Model of the plasma collector and source sheath in potential
Numerical solution for µ = 5 ∗ 10−4 gives us: τ 10 100 1000
α 1,25 1,81 2,38
ψp -0,09 -0,21 -0,25
ψc -0,62 -0,69 -0,74
Table 1: Numerical results
12
3.3
Computer simulation result
When we simulate the same problem in xpdp1 we get the following results:
(a) Potential
(b) Particle density
Figure 5: Resulting potential and particle densities from simulation We see that theoretical model and plasma simulation for bounded plasma agree. The last step in the research is to measure plasma potential in laboratory and if data suits above mentioned methods, we can be pretty sure we know what is going on in plasma in our system.
13
In lower figure 6 is particle distribution through our system depending on their velocities. In figure 6(a) and 6(b) are shown good simulation and in figure 6(c) and 6(d) are results of inappropriate chosen parameters - if we choose too large time step or not enough particle cells for given system length (or possibly both).
(a) Electron velocity
(b) Ion velocity
(c) Electron velocity with bad param- (d) Ion velocity with bad parameters eters
Figure 6: Resulting particle distribution regarding velocity along x-axis for ions and electrons. In upper case is the real physical phenomena, and the lower graphs result from numerical effect.
14
4
Conclusion
I introduced basic one dimensional computer simulation in plasma physics. We saw, how particle flux equilibrium is achieved, although due to lack of collisions, we don’t have any thermodynamical relaxation (besides artificial re-Maxwellization through reflux). At this time more realistic simulations are made and constantly upgraded. Simulation nowadays are not limited to one dimension, interaction between particles is not only through electromagnetic field but also by direct collisions. Various cross sections for collisions are determined and integrated into simulations. It is also taken into account collisions of charged particles with neutrals, and ionization within plasma. Next step is to consider secondary emission after particle collide into electrode and knocks out an electron. To shorten simulation time and improve accuracy variable grids are introduced - with fewer points at constant potential and more points at stirred up potential. To speed up simulation and capability to simulate more particles, programs are adjusted to run simultaneous on several computer processors (clusters). New numerical methods can also improve the simulation performance.
15
References [1] R.W.Hockney, J.W.Eastwood, Computer simulation using particles (IOP Publishing Ltd, 1999) [2] C.K.Birdsall, A.B.Langdon, Plasma physics via computer simulation (McGrrawHill, Inc., 1985) [3] http://www.geocities.com/letapk/pic.html , april 2008 [4] http://ptsg.eecs.berkeley.edu/ , april 2008 [5] P.C.Stangeby, The Plasma Boundary of Magnetic Fusion Devices (Taylor & Francis, 2000) [6] L.A.Schwager, C.K.Birdsall, Collector and source sheats of finite ion temperature plasma (Physics of Fluids B: Plasma Physics, 1990) [7] O.Buneman, Dissipation of Currents in ionized media (Phys. Rev. 115:503, 1959) [8] C.K.Birdsall, W.B.Bridges, Space-charge instabilities in electron diodes and plasma converters (J. Applied Physics, vol. 32, no. 12, pp. 2611-2618, Dec. 1961) [9] C.K.Birdsall, D.Fuss, Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation (J. Comput. Phys. 3, 494–511, 1969)
16
