Computational Plasma Physics Kinetic modelling: Part 1 & 2 W.J. Goedheer

1. Computational Plasma Physics Kinetic modelling: Part 1 & 2 W.J. Goedheer FOM-Instituut voor Plasmafysica Nieuwegein, www.rijnh.nl

2. What are kinetic methods and when do they apply Kinetic methods retain information on the velocity distribution (hydrodynamic/fluid methods first integrate over velocity space) Needed when distribution is non-Maxwellian Kinetic methods are to be preferred when “mfp> L” or “coll>T” mfp and coll depend on densities and cross-sections But what are L and/or T?? Examples from (plasma) physics?

3. Kinetic models: non-Maxwellian Collisions electrons mainly with neutral species Low degree of ionization Effective cooling of parts of the energy distribution function Counteracted by Coulomb collisions at high degree of ionization (>10%) E

4. T P High Pow High T Low Pow Low T L T x t Variations in space and time Boundary layers Transition layers Transient phenomena Switching on Modulation

5. Power modulated discharges Observation in experiments UU) optimum in deposition rate Modulate RF voltage (50MHz) with square wave (1 - 400 kHz) RF

6. Modulated discharge in SiH4 Results from a PIC/MC calculation: Cooling and high energy tail

7. Examples of Ion Energy DF at grounded electrode From Th. Bisschops, Thesis TU/e, 1987 Interaction between E-field and ion motion does not result in a shifted Maxwellian

8. Kinetic models: strong spatial variation Very low pressures: L = size of vessel (applies for [e,i,n]) Space charge boundary layer: L = Debye length (applies for [e,i]) Micro-structures (etched trenches): L = size of structure (applies for [e,i,n]) Shocks: L = extension of shock (applies for [e,i,n]) There may be a difference between momentum loss and energy loss

9. Kinetic models: strong temporal variation Microwave discharges / high frequency RF discharges (applies for [e,(i)]) Start-up of discharges (applies for [e,i]) There may again be a difference between momentum loss and energy loss

10. Methods based on direct solution of the Boltzmann equation Tricks to solve BE: Use symmetry if present Expand f in some small parameter

11. A method especially suitable for electrons Electrons have a low mass  high momentum loss in collisions  energy loss in inelastic collisions Elastic scattering  redistribution over a sphere in velocity space Small deviation from isotropic f in the direction of the average velocity Therefore: expansion in Legendre polynomials Pn(cos ) with  the angle between average and actual velocity f = f0(v) + f1(v)P1(cos ) + f2(v)P2(cos )+……. Note: “Amplitudes” depend on absolute value velocity and vary in space and time

12. f0=v*exp(-v) f1 cos()=0.1v*exp(-v/2)cos() An example: f0+f1cos() f0+f1cos()

13. How to calculate the amplitudes fn Elastic collisons move electrons from outside in Transport from Neighbouring volume V V+dv Shift in velocity due to electric field (eEDt/m) Net change in Dt

14. How to calculate the amplitudes fn Balance of number of particles in shell between v and v+dv in dxdydz Transport in real space Transport in velocity space Effect of collisions

15. How to calculate the amplitudes fn Balance of momentum in shell between v and v+dv in dxdydz Transport in real space Transport in velocity space Effect of collisions Note that f1 is a vector, directed along the average velocity (=0)

16. Cutting off at f1: The Lorentz approximation For elastic collisions with atoms/molecules, with mass M:

17. Special case: homogeneous, steady state Temperature gas is zero Constant electric field, average velocity (f1) along E

18. Special case: homogeneous, steady state Solution: Backward integration, tri-diagonal system …

19. Special solutions: Reduced electric field r=-1 ; s=2 : Maxwell r= 0 ; s=4 : Druyvesteyn Druyvesteyn has less energetic electrons

20. Inelastic collisions • Couple parts of the distribution function that are far apart • Example: Excitation • electron looses excitation energy (a few to >10 eV) • electron is set back in velocity Source proportional to vf0(v)NMexc(v) Same holds for ionization: Energy new electrons to be specified

21. Inelastic collisions: two T distribution Noble gases have high first excitation energy For lower energies only elastic energy losses: slow decay of f with v For higher energies large energy losses: fast decay of f with v Resulting distribution is characterised by two ”temperatures” Eexc Eion

22. Ui-du Ui+du Eion+2(Ui-du) Eion+2(Ui+du) Inelastic collisions: ionization In ionization Eion is lost Suppose remaining energy equally divided How many electrons arrive between Ui-du and Ui+du • Ui-du < (U-Eion)/2 < Ui+du • Eion+2Ui-2du <U <Eion+2Ui+2du • So factor 2 from energy range + factor 2 from new electron: • 4f0(u)u1/2ion(u) In steady state problems: new electrons neglected, Usually this has only a minor influence

23. An example, SiH4/H2, with inelastic collisions EEDFs with 4eV av. Energy in SiH4/H2: non-Maxwellian

24. Some quantities (assuming f0 normalized)

25. Use of this approach in modelling Local field approximation: Everything expressed in local E/N-field mobility and diffusion coefficients reaction rates (ionization, excitation) average energy Mean energy approximation: Use solution for various E/N-fields to construct table: (mobility, diffusion coefficient, rates) all as a function of the average energy (cf. table as function of temperature for Maxwellian f) Use fluid energy balance to obtain av.energy in simulation

26. Use of the mean energy approximation Homogeneous gas of given composition, Nb1...bn EEDF from Boltz.Eqn. Homogeneous electric field, constant in time Mobility (e), Diffusion (De) EEDF Average energy (  1.5 kTe) Reaction rates for processes Kproc (E,Nb1,Nb2,..Nbn) Combine results in table for Kproc () , e(), De()

27. Modelling the electrons + Look-up table

28. One step further: time dependent E-field Important characteristic times: Loss of momentum: goes very fast f1 is in equilibrium with E-field Loss of energy: Only fast in case of inelastic collision f0 can be out of equilibrium Example: reaction of f0 in SiH4/H2 : E0cos(t) : behavior depends on ratio  and loss frequencies

29. High frequency: smaller excursion f0 Collision frequency  pressure Therefore: normalized to 1 Torr 10% SiH4, 90% H2 E=Emcos(t), f0 at Em, Em/2, 0, -Em/2 (1,2,3,4) Energy loss Momentum loss From Capitelli et al.: Pl. Chem. Plasma Proc. 8 (1988) 399-424

30. Time dependent, spatially inhomogeneous E field Is possible in principle, but: More than 1 spatial dimension would take too much CPU time Really steep gradients (sheaths) require fn with n>1 Solution: Monte Carlo methods Account in principle for all effects

31. Example: v*f0 in Nitrogen E=3.6*104(x/L)5(10.8sin(t)), -L<x<L =2*80 MHz

32. Electron-electron collisions Electrons efficiently share energy in elastic collisions  Collisions try to establish Maxwell distribution More sophisticated operators conserve momentum and energy

33. 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 Examples: Electrons in a homogeneous electric field Follow sufficient electrons for a sufficient time Obtain distribution over velocities etc.  f0,f1 Positive ions in plasma boundary layer (ions have trouble loosing momentum)

34. Monte Carlo methods: Equation of motion Leap-frog scheme Alternative: Verlet scheme

35. 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

36. Monte Carlo methods: Boris for B-field Equivalent scheme (J.P.Boris), (proof: substitution):

37. 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) (Rn is random number,0<Rn<1) Number of collisions: NMtot v= 1/ per second. Time to next collision: Tcoll=-* ln(1-Rn)

38. 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 When there is a collision: Determine which one: new random number

39. 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

40. ’s normalized to maximum: Draw random number Max 1+2+3+..N+ 0 v*tot v*0 1+2+3+..N v*3 v* 1+2+3 v*2 1+2 1 v*1 v  Monte Carlo methods: Null Collision

41. 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

42. Monte Carlo methods: Effect of collision v1,v2 velocities in lab-frame prior to collision, w1,w2 in center of mass system

43. 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:

44. 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

45. Monte Carlo methods: Applicability • Examples where MC models can be used are: • motion of electrons in a given electric field in a gas (mixture), see practicum • motion of positive ions through a RF sheath (given E(r,t))

46. Monte Carlo methods: Applicability • Main deficiency of Monte Carlo: 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

47. 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”

48. Particle-In-Cell plus Monte Carlo: Bi-linear interpolation xs zi+1=(i+1)z xs, qs=eNs zs zi=iz xi=ix xi+1=(i+1)x i:=i+(xi+1-xs)qs/x i+1:=i+1+(xs-xi)qs/x xj=jx xj+1=(j+1)x ij:=ij+(zi+1-zs) (xj+1-xs) qs/(x z)

49. Particle-In-Cell plus Monte Carlo: Solution of Poisson equation 2 Boundary conditions on electrodes, symmetry, etc. Electric field needed for acceleration of particle: (bi)linear interpolation, field known in between grid points

50. Move particles F v  x Check loss at the walls Interpolate field to particle Collision new v Solve Poisson equation Interpolate charge to grid Particle-In-Cell plus Monte Carlo: Full cycle, one time step