Introduction to the DSMC / PIC-MC software package

The simulation of transport phenomena is an important method to understand the dynamics in low-pressure deposition reactors. This includes the flow of the process gas as well as the transport of reactive gases and precursors. The most general transport equation in classical mechanics is the Boltzmann transport equation

$$\frac{\partial f}{\partial t} + \vec{v}\,\nabla_x f + \frac{\vec{F}}{m}\,\nabla_v f(\vec{x},\vec{v},t) = \left.\frac{\partial f}{\partial t}\right|_{\mathrm{collision}}$$

It describes the time evolution of a particle distribution function $f(\vec{x},\vec{v},t)$ in a 6-dimensional phase space (space and velocity). Such high-dimensional spaces are usually difficult to handle by numerical simulation, e.g. by the finite element method, because the number of mesh elements would increase with the 6th power of the geometric scale. A perturbation analysis (see e.g. Chapman and Cowling (1970)) shows that for a low value of the mean free path (see section Numerical constraints) in comparison to characteristic geometric dimensions of the flow, the Boltzmann equation can be approximated by the Navier Stokes equations, which merely treat the density, velocity and pressure as flow variables in the three-dimensional space. For the Navier Stokes equations, a large variety of solvers based on Finite Element or Finite Volume Method (FEM, FVM) are available as commercial or open-source codes. A prominent FVM based solver for continuum mechanics is OpenFOAM, whereof different versions are maintained either by the OpenFOAM Foundation or by the ESI Group. Both versions differ slightly with respect to the included solvers and auxiliary libraries.

The validity of the Navier Stokes equations depends on the mean free path $\lambda$ with respect to characteristic geometric dimensions $d$, which can be expressed as the so-called Knudsen number

$$\mathrm{Kn} = \frac{\lambda}{d}$$

As shown in the following figure, continuum approaches are only valid for Knudsen numbers well below 0.1. For higher Knudsen values, kinetic simulation codes must be used instead. In case of very large Knudsen numbers $\mathrm{Kn} \gg 1$, the flow can be considered as collisionless; in this case it is possible to use a ray-tracing approach for solving the Boltzmann equation via a particle based approach. An example of a software which solves particle transport under molecular flow conditions is MolFlow developed at CERN. A more detailed discussion on thin film deposition methods and simulation codes is given in Pflug et al. (2016).

Figure 1.1: Kinetic and continuum solvers for transport phenomena and their validity with respect to the flow characteristic (adapted from Bird (1994))

For an intermediate range $0.1 < \mathrm{Kn} < 100$ of the Knudsen number, neither a continuum approach nor a solver based on collisionless flow is viable. With respect to thin film deposition technology, most Physical Vapor Deposition (PVD) and some of the Chemical Vapor Deposition (CVD) processes fall into this intermediate range with respect to their process conditions. Thus, their proper modeling requires to use a kinetic solver of the Boltzmann equation including collision term.

Other than a fully-fledged FEM or FVM based solution of the Boltzmann equation, a statistical approach called Direct Simulation Monte Carlo (DSMC) - method (see Bird (1994)) is a solution method with feasible computational effort. It is based on following simplifications:

• All flow species are represented by super particles, a super particle has the same mass and properties as real gas molecules of the respective species but represents a number NREAL of real particles. The scaling factor is necessary as otherwise the number of particles to be simulated would be too high to fit into usual high performance computing hardware.
• For collisions between particles, only binary collisions are considered.
• For performance reasons, the particle movement and collision are treated not simultaneously but in two separate routines, which are invoked subsequently for each time cycle.

The 3rd simplification requires that the product $\delta t \times v$ of time step and mean thermal velocity $v$ of all treated species is below the cell spacing of the simulation grid. A typical schedule of the time step in a DSMC solver can be seen in the following figure by following the blue boxes.

Figure 1.2: Scheme of a time cycle in the Direct Simulation Monte Carlo (DSMC) and Particle-in-Cell Monte Carlo (PIC-MC) methods

The DSMC method only treats the transport of neutral particles and hence is suitable to model rarefied gas flows, evaporation or thermal CVD processes. If the transport of charged particles becomes relevant as in plasma processes, the electric and magnetic fields have to be additionally considered. The approximation used in this software is referred to as Vlasov Poisson equation, which means that the force term in the equation above is replaced by the Lorentz force,

$$\frac{\partial f}{\partial t} + \vec{v}\,\nabla_x f + \frac{q}{m}\left(\vec{E} + \vec{v}\times\vec{B}\right)\nabla_v f(\vec{x},\vec{v},t) = \left.\frac{\partial f}{\partial t}\right|_{\mathrm{collision}}$$

and the electric potential $\varphi$ is derived from the Poisson equation according to the actual charge density distribution $\rho$ and the relative dielectric permeability $\varepsilon_r$ of the background.

$$\Delta\varphi + \frac{\nabla\varepsilon_r\,\nabla\varphi}{\varepsilon_r} = -\frac{\rho}{\varepsilon_0\,\varepsilon_r}$$

In a self-consistent Particle-in-Cell Monte Carlo (PIC-MC) - simulation, the last two equations have to be fulfilled simultaneously, while in the actual iterative implementation they are solved sequentially

• as shown by the large cycle in Figure 1.2. This imposes additional numerical constraints on the simulation parameters such as cell spacing, time step etc. as elaborated further in Numerical constraints. A detailed discussion of the PIC-MC method is given in the monography of Birdsall and Langdon (2005).