User Tools

Trace: • num_constraints

Table of Contents

Numerical Constraints

Limitations to the time step width

Cell spacing

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:

$\Delta t \cdot v < \Delta s$.

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:

$\dfrac{m\,v^2}{2\,q} = E_\textrm{kin}[\textrm{eV}] = U_\textrm{max}[\textrm{V}] \quad\Rightarrow\quad v = \sqrt{\dfrac{2\,U\,q}{m}} \left[\dfrac{\textrm{m}}{\textrm{s}}\right]$,

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] vElectron [m/s] vElectron · Δt [μm]
100 5.932 · 106 59.32
200 8.389 · 106 83.89
300 1.027 · 107 102.7
400 1.186 · 107 118.6
500 1.326 · 107 132.6

Collision frequency

According to the Wikipedia1) 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 :

$Z = n_A\,n_B\,\sigma_{AB}\sqrt{\dfrac{8\,k_B\,T}{\mu_{AB}\,\pi}}\,\,[m^{-3}\,s^{-1}]$

where

  • $n_A\,[m^{-3}]$ is the number density of species A
  • $n_B\,[m^{-3}]$ is the number density of species B
  • $\sigma_{AB}\,[m^2]$ is the collision cross section of A and B
  • $k_B\,[m^2\,kg\,s^{-2}\,K^{-1}]$ is the Boltzmann's constant
  • $T\,[K]$ is the temperature
  • $\mu_{AB}\,[kg]$ 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.

Electron collision frequency

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:

$\lambda = N\,\pi\,d^2\sqrt{\dfrac{8\,k_B\,T}{m\,\pi}}$.

Whereas N is the number density, d the effective atom radius, kB the Boltzmann constant, T the gas temperature and m the atom mass. Transfered to the electron collision frequency in a partly ionized plasma one could use the above estimation with:

  • πd2 → σ, the electron-neutral collision cross section with the primary gas species
  • T → Te, the electron temperature whereas 1 eV = 11604.5 K
  • m → me, the electron rest mass of 9.1091 · 10-31 kg

Thus, for a H2-plasma with TH2 = 500 K, Te = 5 eV and an approximated σ = 10-19 m2 one gets:

  • $\lambda=1.0838\cdot 10^{10}\,\textrm{Hz at }P_{H_2} = 500\,\textrm{Pa}$
  • $\lambda=2.1677\cdot 10^{10}\,\textrm{Hz at }P_{H_2} = 1000\,\textrm{Pa}$
  • $\lambda=4.3355\cdot 10^{10}\,\textrm{Hz at }P_{H_2} = 2000\,\textrm{Pa}$

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

Plasma frequency

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

$\omega_p = \sqrt{\dfrac{n\,e^2}{\epsilon_0\,m_\textrm{e}}}$

  • n = Electron density
  • e = 1.6021 · 10-19 As
  • ε0 = 8.8542 · 10-12 F/m
  • me = 9.1091 · 10-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]
1015 1.78 · 109 5.605 · 10-10
1017 1.78 · 1010 5.605 · 10-11
1019 1.78 · 1011 5.605 · 10-12

Debye length

The Debye length LD 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 LD. For electrons, the Debye length is given by:

$L_D = \sqrt\frac{\epsilon_0 k_B T_e}{n_e e^2}} \simeq 69.01\cdot\sqrt{\frac{T_e}{n_e}} \simeq 6070\cdot\sqrt{\frac{E\,[\textrm{eV}]}{n_e}}$.

Here, ε0 is the vacuum permittivity, kB the Boltzmann constant, Te the electron temperature, ne the electron density and e the electron charge.

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 π LD. In [Tskhakaya2007]3) a constraint of:

$\delta_x \le 3.4 L_D$

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 1017/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 1015/m³ is what we often have in simulation (especially 3D) as feasible compromise between realistic power density and computation time.

As an example, for ne = 1015m-3 and Te=5 eV the constraint 3.4 LD is computed to 1.459 mm.

If the simulation output shows “Failure” at the Debye length check this is merely an informative message based on the assumed plasma density of 1015/m³ or 1017/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.

Gyrokinetics

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:

$R_{\scriptscriptstyle\textrm{Larmor}}=\dfrac{m_\textrm{e}\,\left|\vec{v}\right|}{q_\textrm{e}\,\left|\vec{B}\right|}$, with

  • me = 9.1091 · 10-31 kg
  • qe = 1.6021 · 10-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 · 106 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.

Larmor-Rotation at Different Time-Step Widths
Timestep Width dt=1e-11 Timestep Width dt=10e-11
Time Step Widths = 10-11 [s] Time Step Widths = 10-10 [s]

Limitations to the cell spacing

Isotropic vs. anisotropic cell spacing

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:

$Q_{\scriptscriptstyle\textrm{Node}}(t) = f_{\scriptscriptstyle\textrm{Charge}}(x,t)\,f_{\scriptscriptstyle\textrm{Node}}(x) = \int{\Pi_{\Delta x}(x-v_x\,t)\,\Pi_{\Delta x}(x)\,dx}$.

The resulting plot for this mathematical operation can be seen below. Now if the particle's speed vx is increased or the cell width Δx is decreased the variation of Q(x,t) over time increases:

$\left.\dfrac{Q(t)}{dt}\right|_{(v_x,\,\Delta x)} < \left.\dfrac{Q(t)}{dt}\right|_{(2\,v_x,\,\Delta x)}\quad\textrm{and}\quad\left.\dfrac{Q(t)}{dt}\right|_{(v_x,\,2\,\Delta x)} < \left.\dfrac{Q(t)}{dt}\right|_{(v_x,\,\Delta x)}$.

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

$\dfrac{Q(t)}{dt}\approx \dfrac{E(t)}{dt}$.

1)
https://en.wikipedia.org/wiki/Collision_frequency
2)
LANGDON, A. B.: Effects of the Spatial Grid in Simulation Plasmas. In: Journal of Computational Physics 6 (1970), p. 247-267
3)
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

Page Tools

SIMULATION

Site Tools

Imprint |  Data Protection Policy