One limitation for the maximum time step width is given due to the implemented Monte Carlo algorithm. According to their velocities particles must stay within a grid cell for at least on time step. Otherwise the collision probabilities are not evaluated correctly, i.e. a »cell hopping« occurs. Hence, the product of time step and particle velocity must not exeed the cell spacing:

.

The velocity of charged particles can be estimated from the maximum voltage in the simulation setup. If neglecting the energy loss due to inelastic collisions, the voltage potential is applied as kinetic energy to the charged particle:

,

with *q* = 1.60211 · 10^{-19} [As]. The table below lists electron velocities for various energies and the accordant distance covered at a time step width of **10 ^{-11} [s]**:

Energy [eV] | v_{Electron} [m/s] | v_{Electron} · Δt [μm] |
---|---|---|

100 | 5.932 · 10^{6} | 59.32 |

200 | 8.389 · 10^{6} | 83.89 |

300 | 1.027 · 10^{7} | 102.7 |

400 | 1.186 · 10^{7} | 118.6 |

500 | 1.326 · 10^{7} | 132.6 |

According to the Wikipedia^{1)} the collision frequency Z can be descibed as the rate of collisions between two species in an given volume and time. Regarding both species A and B as hard spheres, Z is expressed as :

where

- is the number density of species A
- is the number density of species B
- is the collision cross section of A and B
- is the Boltzmann's constant
- is the temperature
- is the reduced mass of A and B

Since we are interested in the lower bound of the time step, only the highest cross section is of interest which is usually the reaction of the main gas species with itself.

According to http://www.uccs.edu/~tchriste/courses/PHYS549/549lectures/plasma.html the collision frequency for neutrals in an ideal gas can be estimated from:

.

Whereas *N* is the number density, *d* the effective atom radius, *k _{B}* the Boltzmann constant,

*πd*, the electron-neutral collision cross section with the primary gas species^{2}→ σ*T → T*, the electron temperature whereas 1 eV = 11604.5 K_{e}*m → m*, the electron rest mass of 9.1091 · 10_{e}^{-31}kg

Thus, for a H_{2}-plasma with *T _{H2}* = 500 K,

For a valid simulation, time step width must be well below the reciprocal of the collision frequency.

The plasma frequency *ω _{p}* is the angular frequency of the electron oscillations around the ions. It depends only on electron mass and density:

*n*= Electron density*e*= 1.6021 · 10^{-19}As*ε*= 8.8542 · 10_{0}^{-12}F/m*m*= 9.1091 · 10_{e}^{-31}kg

In order to properly resolve the plasma oscillations, the electron time step should be well below the respective oscillation period e.g., dt < 0.1 *ω _{p}*, as can be estimated from the following table:

Electron density [m^{-3}] | Plasma frequency [s^{-1}] | Oscillation period [s] |
---|---|---|

10^{15} | 1.78 · 10^{9} | 5.605 · 10^{-10} |

10^{17} | 1.78 · 10^{10} | 5.605 · 10^{-11} |

10^{19} | 1.78 · 10^{11} | 5.605 · 10^{-12} |

The Debye length *L _{D}* describes the electrical shielding capability of the plasma. The potential difference induced by a local excess charge will drop to 1/e of the original potential at the distance of

.

Here, *ε _{0}* is the vacuum permittivity,

In order to properly resolve the electric potential distribution in the simulation grid, the cell size has to in the same order of magnitude as the Debye length. From [Langdon70]^{2)} follows, that the cell spacing shall be less than *π L _{D}*. In [Tskhakaya2007]

is suggested. This is what we also adopted in our constraint check message for the Debye length on plasma simulation startup. If this constraint is broken, artificial numerical oscillations will take over leading to “artificial heating”, i.e. more energy will be gained in the plasma than actually put into the system via the electrodes, and the simulation will become unstable.

Our simulation tool cannot predict the plasma density of any given simulation case. Therefore, we give some estimations for typical cases:

- We assume 5 eV for electrons as this is a very typical electron temperature in magnetron sputtering discharges.
- A plasma density of 10
^{17}/m³ corresponds to an experimentally realistic density (still with moderate power density, e.g. 50-100 W for a 3 inch cathode). - A plasma density of less than 10
^{15}/m³ is what we often have in simulation (especially 3D) as feasible compromise between realistic power density and computation time.

As an example, for *n _{e}* = 10

If the simulation output shows “Failure” at the Debye length check this is merely an informative message based on the assumed plasma density of 10^{15}/m³ or 10^{17}/m³ and electron temperature of 5 eV. If the actual simulation setup is different from that it may be still numerically valid. If in doubt it helps to check if two simulation runs with slightly changed cell spacing yield the same overall results.

The implemented leap-frog algorithm from Boris is very robust at comparatively small computation costs. Nevertheless if the chosen time step witdh is too large, gyrokinetic trajectories could not be resolved sufficiently as seen in the pictures below. For a proper estimation use the following equation:

, with

*m*= 9.1091 · 10_{e}^{-31}kg*q*= 1.6021 · 10_{e}^{-19}As

This formular yields the path radius of an electron at constant velocity *v* perpendicular to a homogenous magnetic field *B*. **We still lack informations on how the time step related granularity is influencing the
plasma simulation and the results respectively.** However, for a proper simulation of gyrokinetic plasmas it should
not exeed a certain limit. Apparant from the pictures below, at magnetic field strengths around 100 mT and electron
velocities around 5 · 10^{6} m/s, the time step width for electrons should not exeed 10^{-10} s. These values are commonly valid for magnetron discharges in PVD sputtering processes on which the PICMC simulation environment is focused on.

Within the Particle-in-Cell (PIC) algorithm charged particles are mapped on the cell nodes by a weighting scheme. Then, in a successive computation step the Poisson Equation of the electric field can be solved consistently to the weighted charge distribution. The PIC weighting scheme in general delocalized a point charge to a finite sized space charge with an appropriate shape function. The charge on a certain node is then given by the fraction of the space charge within the cell volume.

In a homogenous grid the particle shape is equivalent to the cell spacing at every point. Thus, for a bilinear (trilinear) weighting scheme, the charge on a certain cell node can be described as the convolution of two square pulses over time:

.

The resulting plot for this mathematical operation can be seen below. Now if the particle's speed *v _{x}* is increased

.

If the charge on a nodes varies over time so does the electric field resulting from the solution of the poisson equation:

.

LANGDON, A. B.: Effects of the Spatial Grid in Simulation Plasmas. In: Journal of Computational Physics 6 (1970), p. 247-267

TSKHAKAYA, D.; MATYASH, K.; SCHNEIDER, R. & TACCOGNA, F.: The Particle-In-Cell Method. In: Contributions to Plasma Physics 47 (2007), No. 8-9, p. 563-594

**How-to**

**Modeling 101**

**Examples**

**Expert corner**

**Software internals**