Transcript Particle-In-Cell plus Monte Carlo

Computational Plasma Physics
Kinetic modelling: Part 2
W.J. Goedheer
FOM-Instituut voor Plasmafysica
Nieuwegein, www.rijnh.nl
Monte Carlo methods
Principle: Follow particles by
- solving Newton’s equation of motion
- including the effect of collisions
- collision: an event that instantaneously changes the velocity
Note: The details of a collision are not modeled
Only the differential cross section + effect on energy is used
Example: Electrons in a homogeneous electric field
Follow sufficient electrons for a sufficient time
Obtain distribution over velocities etc.  f0,f1
Monte Carlo methods: Equation of motion
Leap-frog scheme
Monte Carlo methods: B-field
Problem with Lorentz force: contains velocity, needed at time t
Solution: take average
The new velocity at the right hand side can be eliminated by taking the
cross product of the equation with the vector
Monte Carlo methods: Boris for B-field
Equivalent scheme (J.P.Boris), (proof: substitution):
Monte Carlo methods: Collisions
Number of collisions: NMtot = 1/ per meter.
 (x) = (0)*exp(- NMx) = (0)*exp(-x/)
dP(x)=fraction colliding in (x,x+dx)=exp(-x/)(1-exp(-dx/))=(dx/)exp(-x/)
 P(x)=(1-exp (-x/))
Distance to next collision: Lcoll=-*ln(1-Rn)
Number of collisions: NMtot v= 1/ per second.
Time to next collision: Tcoll=-* ln(1-Rn)
(Rn is random number,0<Rn<1)
Monte Carlo methods: Collisions
Another approach is to work with the chance
to have a collision on vt: Pc=vt/
Ensure that vt<< to have no more than one collision per timestep
 Effect of collision just after advancing position or velocity
 introduces only small error
1.0
no collision
colliding once
colliding twice
colliding> twice
0.9
0.8
When there is a collision:
fraction
Determine which one: new random number
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.5
1.0
1.5
L/
2.0
2.5
3.0
Monte Carlo methods: Null Collision
Problem: Mean free path is function of velocity
Velocity changes over one mean free path
Solution: Add so-called null-collision to make v*tot independent of v
Null-collision does nothing with velocity
Mean free path thus based on Max (v*tot)
Is rather time-consuming when v*tot peaks strongly
Monte Carlo methods: Null Collision
Max
v*tot
v*0
v*
v*3
v*2
v*1
v
’s normalized to maximum:
Draw random number
1+2+3+..N+ 0
1+2+3+..N
1+2+3
1+2
1
Monte Carlo methods: Effect of collision
Determine effect on velocity vector
Retain velocity of centre of gravity
Select by random numbers two angles of rotation for relative velocity
Subtract energy loss from relative energy
Redistribute relative velocity over collision partners
Add velocity centre of gravity
Monte Carlo methods: Effect of collision
v1,v2 velocities in lab-frame prior to collision, w1,w2 in center of mass system
Monte Carlo methods: Effect of collision
A collision changes the size of the relative velocity if it is inelastic
A collision rotates the relative velocity
Two angles of rotation:    and   
 usually has an isotropic distribution: =Rn*
 has a non-isotropic distribution
Hard spheres:
Monte Carlo methods: Rotating the relative velocity
Step 1: construct a base of three unit vectors:
Step 2: draw the two angles
Step 3: construct new relative velocity
Step 4: construct new velocities in center of mass frame
Step 5: add center of mass velocity
Monte Carlo methods: Applicability
Examples where MC models can be used are:
- motion of electrons in a given electric field in a gas (mixture)
- motion of positive ions through a RF sheath (given E(r,t))
Main deficiency: not selfconsistent
- electric field depends on generated net electric charge distribution
- current density depends on average velocities
- following all electrons/ions is impossible
Way out: Particle-In-Cell plus Monte Carlo approach
Particle-In-Cell plus Monte Carlo: the basics
-Interactions between particle and background gas are dealt with only in collisions
-this means that PIC/MC is not! Molecular Dynamics
-each particle followed in MC represents many others: superparticle
-Note: each “superparticle” behaves as a single electron/ion
-Electric fields/currents are computed from the superparticle densities/velocities
-But: charge density is interpolated to a grid, so no “delta functions”
Particle-In-Cell plus Monte Carlo:
Bi-linear interpolation
xs, qs=eNs
xi=ix
xi+1=(i+1)x
i:=i+(xi+1-xs)qs/x
i+1:=i+1+(xs-xi)qs/x
zi+1=(i+1)z
xs
zs
zi=iz
xj=jx
xj+1=(j+1)x
ij:=ij+(zi+1-zs) (xj+1-xs) qs/(x z)
Particle-In-Cell plus Monte Carlo:
Solution of Poisson equation
Boundary conditions on electrodes, symmetry, etc.
Electric field needed for acceleration of particle:
(bi)linear interpolation, field known in between grid points
Particle-In-Cell plus Monte Carlo:
Full cycle, one time step
Move particles
F v  x
Interpolate field
to particle
Solve Poisson
equation
Check loss
at the walls
Collision
new v
Interpolate
charge to grid
Particle-In-Cell plus Monte Carlo:
Problems
Main source of problems: Statistical fluctuations
Fluctuations in charge distribution: fluctuations in E
average is zero but average E2 is not  numerical heating
Sheath regions contains only few electrons
Tail of energy distribution contains only few electrons
large fluctuations in ionization rate can occur
Particle-In-Cell plus Monte Carlo:
Problems
Solutions:
-Take more particles (NB error as N-1/2 ) , parallel processing!
-Average over a long time
-Split superparticles in smaller particles when needed
requires a lot of bookkeeping, different weights!
Particle-In-Cell plus Monte Carlo:
Stability
Plasmas have a natural frequency for charge fluctuations:
The (angular) Plasma Frequency:
And a natural length for shielding of charges:
The Debye Length:
Stability of PIC/MC requires:
Power modulated discharges
Modulate RF voltage (50MHz)
with square wave (1 - 400 kHz)
Observation in experiments UU)
optimum in deposition rate
Modulated discharges
Results from a PIC/MC calculation: Cooling and high energy tail
1-D Particle-In-Cell plus Monte Carlo Simulation
of a dusty argon plasma
Dust particles with a homogeneous density distribution are present in two layers
This resembles certain experiments done under micro-gravity
Dust particles do not move, they only collect and scatter plasma ions and electrons
The charge of the dust results from the collection process
The charge of the dust is defined on the grid needed for the Poisson equation
1-D Particle-In-Cell plus Monte Carlo Simulation
of a dusty argon plasma
Crystal (21010 m-3)
7.5 m radius
Void
Capture cross section
Scattering:
Coulomb, truncated at d
L/8 L/4
RF
w is energy electron/ion
Charging of the dust upon capture of ion/electron
The total charge is monitored on the gridpoints
Charge of captured superparticle is added to nearest gridpoints
Division according to linear interpolation
Superparticle is removed
Local dustparticle charge is total charge
divided by nr. of dust particles
This number is: density*dz*a2, with a the electrode radius
For Monte Carlo the maximum v is computed for all
available dust particle charges
Null-collision is used
Simulation for Argon, 50MHz, 100mTorr, 70V, L=3cm
dustfree
with dust
1.50x1016
6x1015
N
1.25x1016
N
N
Density (m -3)
4x1015
7.50x1015
5.00x1015
2.50x1015
d
e
+
3x1015
2x1015
1x1015
0.00
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0
0.0
1.0
0.1
0.2
0.3
0.4
x/L
108
0.5
0.6
0.7
0.8
0.9
1.0
x/L
108
107
104
106
EEDF (arb.un.)
105
103
104
103
102
101
101
5
10
15
20
Energy (eV)
25
30
35
40
L/2
L/4
3L/16
L/8
L/16
105
102
0
Vd6V
107
L/2
L/4
L/8
L/16
L/32
106
EEDF (arb.un.)
Density (m -3)
d
N
+
1.00x1016
100
N Q /e
5x1015
e
100
0
5
10
15
20
Energy (eV)
25
30
35
40
Simulation for Argon, 50MHz, 100mTorr, 70V, L=3cm
dustfree
with dust
3.0
3.0
Av. El. Energy
Ion.Rate
Av. El. Energy
2.0
1.5
e
1.0
e
3kT /2 (eV), Ion.Rate (arb.un.)
3kT /2 (eV), Ion.Rate (arb.un.)
Ion.Rate
2.5
0.5
0.0
0.000
0.125
0.250
0.375
0.500
0.625
0.750
0.875
1.5
1.0
0.5
0.125
0.250
0.375
106
0.750
106
105
104
103
0.875
1.000
0
L/16
L/8
3L/16
L/4
L/2
107
IEDF (arb.un.)
107
0.625
108
L/2
L/8
L/16
L/32
0
108
0.500
x/L
109
109
IEDF (arb.un.)
2.0
0.0
0.000
1.000
x/L
105
104
103
102
102
101
101
100
2.5
100
0
5
10
15
20
Energy (eV)
25
30
35
40
0
5
10
15
20
Energy (eV)
25
30
35
40
Simulation for Argon, 50MHz, 100mTorr, 70V, L=3cm
Generation of internal space charge layers
3x1014
Net charge / e
2x1014
An internal sheath is formed
inside the crystal
1x1014
0
-1x1014
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
x/L
20000
Ions are accelerated before
they enter the crystal
15000
Av. Electric Field (V/m)
10000
This has consequences for
the charging + shielding
5000
0
-5000
-10000
-15000
-20000
0.000
0.125
0.250
0.375
0.500
x/L
0.625
0.750
0.875
1.000
Particle-In-Cell plus Monte Carlo:
What if superparticles collide?
Example: recombination between positive and negative ions
Procedure: number of recombinations in t: N+N-Krec t
corresponds to removal of corresponding superparticles
randomly remove negative ion and nearest positive ion
but: be careful if distribution is not homogeneous
A more sophisticated approach: Direct Simulation Monte Carlo
DSMC: Basics
Divide the geometry in cells
Each cell should contain enough testparticles (typically 25)
Newton’s equation: as before, but keep track of cell number
Collisions: choose pairs (in same cell!) and make them collide
Essential: the velocity distribution function is sum of -functions
Only small fraction of pairs collides in one time step
DSMC: Choosing the pairs
Add null collision
Chance of collision of particle i with j is Pc=(Npp/Vcell)*Max(v)t
Number of colliding pairs: n(n-1)* Pc/2
Select randomly particle pairs (make sure no double selection)
See if there is no null collision (again with random number)
Perform the collision
DSMC: An example
Relaxation of a mono-energetic distribution to equilibrium
20000 particles, hard sphere collisions. All particles are in
the same cell. Distribution at various time steps
3200
25
50
75
100
2800
# particles
2400
2000
1600
1200
800
400
0
0
16
32
Energy (arb.un.)
48
64