2.1.3.1 Numerical example

We assume that a turbomolecular pump with a nominal pumping speed of \(S_p = 1.0\) m³/s for Argon at 300 K is attached to a circular surface with a diameter of 133 mm. From that, the following steps are necessary to compute the effective pump factor \(f\):

• The area \(A\) of the flange is \(\pi r^2\) with r=66.5 mm; thus, A = 0.01389 m²

• From to Eqn. ((eq:mean_thermal_speed?)) follows that \(\overline{c} = 399\) m/s for Argon at 300 K

• According to Eqn. ((eq:pumping_speed_orifice?)) the maximum possible pumping speed is 1.386 m³/s

• From that, Eqn. ((eq:pump_factor?)) yields \(f = 72.16\)%

2.1.4.1 Numerical example

In a chamber with the turbomolecular pump described in Section (sec:num_example_sp?), Argon is injected with a flow rate of 500 sccm, and the mean gas temperature is 300 K. According to Eqn. ((eq:pressure_flow_pumping?)), the equilibrium pressure will be

\[p = \Phi\frac{F_{in}}{S_p} = 0.001854\textrm{ (Nm/(sccm s))} \frac{500\textrm{ (sccm)}}{1\textrm{ (m}^3\textrm{/s)}} = 0.927\textrm{ (Pa)}\]

2.2.3.1 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 \times v < \Delta s\]{#eq:constraint_cellspacing}

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:

\[\frac{m\,v^2}{2\,q} = E_\textrm{kin} (\textrm{eV}) = U_\textrm{max} (\textrm{V}) \Rightarrow v = \sqrt{\frac{2\,U\,q}{m}} \left[\dfrac{\textrm{m}}{\textrm{s}}\right]\]{#eq:constraint_charge_velocity}

with \(q = 1.60211 \times 10^-19\) (As). The table below lists electron velocities for various energies and the accordant distance covered at a time step width of \(\delta t=10^{-11}\) (s):

2.2.3.2 Collision frequency

According to the Wikipedia¹ the collision frequency \(nu\) 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, \(nu\) is expressed as

\[\nu = n_A n_B \sigma_{AB} \sqrt{\frac{8 k_B T}{\mu_{AB}\pi}}\,(m^{-3} s^{-1})\]{#eq:collision_frequency}

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.

2.2.3.3 Electron collision frequency

The collision frequency for neutrals in an ideal gas can be estimated from

\[\nu = N\,\pi d^2\sqrt{\frac{8 k_B T}{m \pi}}\]{#eq:collision_frequency_2}

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

• \(\pi d^2 \rightarrow \sigma\) = the electron-neutral collision cross section with the primary gas species
• \(T \rightarrow T_e\) = the electron temperature whereas 1 eV = 11604.5 K
• \(m \rightarrow m_e\) = the electron rest mass of \(9.1091 \times 10^{-31}\) kg

Thus, for a H₂-plasma with T_H2 = 500 K, T_e = 5 eV and an approximated \(\sigma\) = 10^-19 m² one gets:

• \(\lambda=1.0838\times 10^{10}\,\textrm{Hz at }P_{H_2} = 500\,\textrm{Pa}\)
• \(\lambda=2.1677\times 10^{10}\,\textrm{Hz at }P_{H_2} = 1000\,\textrm{Pa}\)
• \(\lambda=4.3355\times 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.

2.2.3.4 Plasma frequency

The plasma frequency \(\omega_p\) is the angular frequency of the electron oscillations around the ions. It depends only on electron mass and density:

\[\omega_p = \sqrt{\frac{n e^2}{\epsilon_0 m_e}}\]{#eq:omega_p}

• \(n\) = Electron density
• \(e\) = Electron charge (\(1.6021\times 10^{-19}\) As)
• \(\epsilon_0\) = Dielectric permittivity (\(8.8542\times 10^{-12}\) F/m)
• \(m_e\) = Electron mass (\(9.1091\times 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., \(\delta t < 0.1 \omega_p\), as can be estimated from the following table:

2.2.3.5 Debye length

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 \(L_D\). 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\times\sqrt{\frac{T_e}{n_e}} \simeq 6070\times\sqrt{\frac{E (\textrm{eV})}{n_e}}\]

Here, \(\epsilon_0\) is the vacuum permittivity, \(k_B\) the Boltzmann constant, \(T_e\) the electron temperature, \(n_e\) 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 (Langdon 1970) follows, that the cell spacing shall be less than \(\pi L_D\). In (Tskhakaya et al. 2007) 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 10¹⁷/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¹⁵/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^15 m^{-3}\) and \(T_e = 5 eV\), the constraint \(3.4 L_D\) 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 10¹⁵/m³ or 10¹⁷/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.

2.2.3.6 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_{\textrm{Larmor}} = \frac{m_e\left|\vec{v}\right|}{q_e\left|\vec{B}\right|}\]{#eq:rlarmor}

with

• \(m_e = 9.1091\times 10^{-31}\) kg
• \(q_e = 1.6021\times 10^{-19}\) As

Eqn. ((eq:rlarmor?)) 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\times 10^6\) m/s, the time step 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 computed at different time steps. {#tbl:larmor_rotation}

Time Step = 10-11 [s]	Time step = 10-10 [s]

2.2.3.7 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_{Node}(t) = f_{Charge}(x,t)f_{Node}(x) = \int dx\, \Pi_{\Delta x}(x-v_x t)\Pi_{\Delta x}\]

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

\[\left.\frac{Q(t)}{dt}\right|_{(v_x, \Delta x)} < \left.\frac{Q(t)}{dt}\right|_{2v_x, \Delta x}\,\,\textrm{and}\,\,\left.\frac{Q(t)}{dt}\right|_{(v_x, 2 \Delta x)} < \left.\frac{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:

\[\frac{Q(t)}{dt}\approx \frac{E(t)}{dt}\]