Table of Contents Author Guidelines Submit a Manuscript
Modelling and Simulation in Engineering
Volume 2011, Article ID 931415, 13 pages
Research Article

Plasma and BIAS Modeling: Self-Consistent Electrostatic Particle-in-Cell with Low-Density Argon Plasma for TiC

1Department of Mathematics, Humboldt-University of Berlin, Unter den Linden 6, 10099 Berlin, Germany
2Department of Physics, Humboldt-University of Berlin, 10099 Berlin, Germany

Received 10 May 2011; Accepted 22 June 2011

Academic Editor: Zeki Ayag

Copyright © 2011 Jürgen Geiser and Sven Blankenburg. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We motivate our study by simulating the particle transport of a thin film deposition process done by PVD (physical vapor deposition) processes. In this paper we present a new model taken into account a self-consistent electrostatic-particle in cell model with low density Argon plasma. The collision model are based of Monte Carlo simulations is discussed for DC sputtering in lower pressure regimes. In order to simulate transport phenomena within sputtering processes realistically, a spatial and temporal knowledge of the plasma density and electrostatic field configuration is needed. Due to relatively low plasma densities, continuum fluid equations are not applicable. We propose instead a Particle-in-cell (PIC) method, which allows the study of plasma behavior by computing the trajectories of finite-size particles under the action of an external and self-consistent electric field defined in a grid of points.

1. Introduction

We motivate our study by simulating a thin-film deposition process that can be done with PVD (physical vapor deposition) processes or sputtering processes.

In the last years, the research in producing high temperature films by PVD (physical vapor deposition) processes in low pressure regimes has increased. Therefore, the interest on standard applications to TiN and TiC are immense but recently also deposition with new material classes known as MAX-phases are important and we are motivated by such applications. The MAX phase is nanolayered terniar metal carbides or nitrids, where M is a transition metal, A is an A-group element (e.g., Al, Ga, In, Si, etc.), and X is C (carbon) or N (nitride).

We present a model for low-temperature and low-pressure plasma. The model is derived as a pathway model, see [1], to achieve the deposition rates of the stoichiometry Ti and C.

We taken into account the drift of the ions in the reactor and the less retardation of molecules, which are not treated by the plasma.

The model is discussed as a mass-conserved transport problem, modeled with convection-diffusion and reaction equations.

Additionally, we have modeled the electric and magnetic field for the particles, which are controlled by the BIAS voltage of the electrostatic field.

We apply the electric and magnetic field in the Lorentz forces to obtain a self-consistent electrostatic particle in cell model. At the target we simulate the deposition rates of the particles Ti and C based on the Monte Carlo simulations, see idea in [2, 3].

We compare with physical experiments, which measure the different ions and molecule rates in the reactor.

The paper is outlined as follows. In Section 2, we present our mathematical model and a possible reduced model for the further approximations. The methods are given in Section 3. Numerical experiments are given in Section 4. In the contents that are given in Section 5, we summarize our results.

2. Mathematical Model

In the following, we discuss the model of the particle transport in the apparatus between substrate and target. In a next subsection, we describe the electromagnetic field that influences the transport of the ionized particles.

2.1. Collision Model: Mean Free Path

The mean free path or average distance between collisions for a gas molecule may be estimated by kinetic theory. If one assumes the gas to be consisted of hard spheres (nonoverlapping spheres), then the effective collision area is given by In time , the area sweeps out the volume , and the number of collisions can be estimated from the number of target molecules () that are in that volume This expression for the mean free path is a good approximation, but it suffers from a significant flaw; it assumes the target objects to be at rest, which is of course physically nonsense. By introducing a relative velocity between the gas objects, whereby the results from the molecular speed distribution of a monoatomic ideal gas (Maxwell Boltzmann distribution). We therefore have the expression The number of molecules per unit volume can be determined from the state equation of the gas If one assumes an ideal gas (noninteraction and nonoverlapping gas particles), one can neglect the so-called higher virial coefficients (). Inserting the state equation for an ideal gas into (4), one gets whereby is the gas constant, and is Avogadro’s number. This is an approximation for mean free path for an atom/molecule of an ideal gas. In our problem, however, we have to calculate the mean free path of an external particle (projectile) which is not a member of the background gas (ideal gas). This can be done by modifying the average relative velocity between projectile and target. This is done in the next part.

2.2. The Mean Relative Velocity between Projectiles and Targets

The background gas is assumed to be Maxwell distributed in velocity (this is motivated by the assumption of an ideal gas). Because of the fact that the background particles being a particle ensemble (with statistically distributed velocities), one can just speak of a mean relative velocity , which can be calculated via where is the three-dimensional Maxwell distribution given by with the abbreviation . A complete derivation of the solution can be found in the Appendix. The result is with (scalar) and . We now want to discuss a few special cases.

If the velocity of the projectile is very small , then , and therefore, the following approximation holds: which gives (2) as expected.

If the targets objects are identically the projectile objects (same mass and same mean velocity), then the following limit holds: which gives the factor and leads to the mean free path of an element of a monoatomic ideal gas (as expected). However, the general expression for the mean free path of a projectile probing into an ideal gas with pressure and temperature is given by There are a few things to say about this expression. First, the main assumption that the background gas (ensemble of target particles) is an ideal gas is just valid within the high-vacuum regime, that is, small target density. Second, the interaction between the projectile and target atoms is assumed to be of hard sphere type, that is, purely geometric interaction. If the projectile is a free particle between the interactions, its Hamilton function reads In this case, one can easily compute . It follows immediately that In appropriate units (atomic units), the scalar reads And therefore, the mean free path in units of cm is given by Eklund [5] used a formula for the mean free path of ions surrounded by an ideal gas of pressure given by In the following Table 1, we present the mean free path for ions at eV and K and gas pressure mbar. We apply an argon gas atmosphere as a target gas in the sputter process.

Table 1: Parameter of the Ions.

In a sputtering process, the ions obey a kinetic-energy distribution as well as an angular distribution at the target. Because of different transport mechanism, the ions loose some extend of their initial kinetic energy. An individual ion within a sputter process can therefore be classified into three groups. First, the ballistic group, which is excelled in the way that any member of the ballistic group travels from the target to the substrate in a straight line because no collisions occur. Second, the transition group is characterized by the observation that the path of the ion is not a straight line, and therefore, the ions of this group undergo some collisions but still retain some of their initial energy. The last group is the thermalized or diffusive group, whereby any member of this group is characterized by a complete loss of its initial kinetic energy. The motion of such an ion is therefore described by a random walk. The typical distances between the target and the substrate are of the order 5–15 cm. Hence, at low argon pressures, we can classify carbon as more or less ballistic, and silicium and titanium as transition or thermalized. One can also see that the formula used by Eklund (2007) [5] is a quite good approximation although it lacks an energy dependency of the mean free path with respect to the ion energy. There are several attempts to achieve an energy dependency in the mean free path. But most of them are more or less physical consistents. For example, Mahieu et al. (2006) [6] use a formula, whereby the energy dependency is arrived by modifying the naive mean free path by multiplying the naive formula with the ion energy. This is of course unphysical because it implies an infinite mean free path at very high ion energies (more precisely, the associate cross-section cannot be normalized, that is, unitarity violation of the cross-section). Our mean free path equation is always finite, and therefore, no violation of unitarity is expected. We hope that our formula for the mean free path will positively be accepted within the community and might help to implement a realistic description of the interactions between particles. In Figure 1, one can see the results from (15) and (16) with respect to the ion energy E (kinetic energy) at an argon pressure of  mbar and a constant temperature of K, whereby the following constants were used, see Table 2.

Table 2: Characteristics: Mass and radius of the elements.
Figure 1: Mean free path of projectiles at argon targets ( mbar and K).
2.3. Electromagnetic Field

For the electromagnetic field, we have the following assumption.

Assumption 1. (i) One decouples between ions and electrons (electrons behave like perfect fluid, due to lower mass).
(ii) One applies Monte Carlo particle-in-cell simulations for the transport of the particles (of the order of 10 mio. particles).
(iii) One applies optimized iterative solver for the Maxwell equations.

We deal with the following equations: the Lorentz force for the particles is given as where is the Lorentz force on each particle, is the electric field (in volts per meter), and is the magnetic flux density (in teslas). Further is the electric charge of the particle (in coulombs), and is the instantaneous velocity of the particle (in metres per second).

We denote × as the vector cross-product.

Further, we have , where is the magnetic permeability and is the magnetic field.

Further, we have where is the scalar field representing the electric potential at a given point. is charge density in space, and is the permittivity of free space (electric constant).

Further, we denote while the ion charge is given in the following equations: with is the charge of ion and density .

The electrons as fluid are given as where is the Boltzmann’s constant, is the mean concentration of charges of the electrons, and is the temperatures of the electrons. Further, is the mean charge of the electrons.

We have the following parameters for the computations of the discretized schemes: (i)debey length: ,(ii)thermal velocity of ions: ,(iii)drift velocity of ions: was varied between 5000 and 9000 m/s, (iv)operator discretization via finite difference scheme (spatial: in units of debey length, temporal: ).

We solve the equation of motion solver with the leap frog method (simplest symplectic integrator).

Further we assume the following:(i)sputter particles are Boltzmann distributed with a mean energy of 2 eV, (ii)angular distribution of sputtered particles is assumed to be Gaussian with mean value of 0 degree and a variance of 10 degrees.

3. Monte Carlo Simulations

In the following, we apply the Monte Carlo simulation, that is based on the collision model for the DC sputtering and the Coloumb model for the HIPIMS sputtering, see [7].

3.1. Angular Distribution

The angular distribution of out-coming particles from the sputter material is modeled by a sine distribution, that is, the relative amount of particles leaving the sputter material perpendicular to its surface () is 1. Differences in the angular distribution between the different species are not modeled but cannot experimentally excluded.

3.2. Ionization Rates and Ion Energy Distribution

The ionization rate of sputtered particles is very low, and thus no influence on the particle distribution is assumed. But in contrast, the particle’s energy seems to be of high importance. Unfortunately, until now no energy distribution for our compound target (TiC) is available. In Figure 2, one can see the ion energy distribution, which is modeled with reference to a Ti-target. One can see that most of the ions are at energies close to 3 eV. In order to simulate the ion transport it is necessary to calculate the velocity of the ions. With, it follows that The energy of the ions is given in units of electron volts (eV), and the mass of the ions is given in atomic units (u). Therefore, one can compute the velocity in units of cm per second by using In two spatial dimensions, one has two velocity components. is the direction angle of the ion (see angular distribution of the ions), and the velocity components can be calculated by

Figure 2: (a) energy distribution and (b) angular distribution of the sputtered species at the target.

Remark 1. In Figure 2, the numerical simulations are shown for a Ti-target. One can see the ion energy distribution in the left figure and can see that most of the ions are at energies close to 3 eV (highest peak). The energy distributions of the sputtered Ti-species are close to very low, such that the possibility to rest on the surface of the deposited material is very high. In Figure 2(b), one can see the angle distribution of the sputtered species, which are highest between 60° and 120°. Such results support the assumptions of a directed transport of the sputtered species to the target. The small peaks at 60° and 120° are negligible and are in the tolerance of the numerical experiments.

Now, we want to apply our two interaction models to DC and high-power pulsed sputtering for TiC. In general, if several independent interaction mechanism can occur, the mean free path is not an additive quantity, but in contrast, the total cross-section is an additive quantity. In order to reduce the computational effort, we decided to use an event-driven Monte Carlo method in contrast to the usually used time-driven Monte Carlo method. It is therefore necessary to determine when the next interaction will occur. If the velocity () and the mean free path () of the particle are known, one can compute the collision frequency by using With the help of the collision frequency, one is able to compute the time interval until the interaction occurs whereby is a random number from a uniform distribution between zero and one. Instead of simulating the trajectory of all particles in a Monte Carlo run with a fixed time step, one can use the above-mentioned formula to adjust the time step. The strategy is as follows: one calculates the time interval for every particle (except the background particles) within a Monte Carlo run (trial) and finds the minimum value. The particle related to the minimum value of will first undergo an interaction. The Monte Carlo time step is set to this minimum value (event-driven MC). After the time step, the specific particle will undergo the interaction, and all other particles just move along their specific trajectory, that is, in the absence of any external forces, the trajectory is just a straight line (this is motivated by the fact that even if external fields are set up, inside the plasma the particles will behave as if they were free, due to the electric conductance of the plasma). If an interaction with the background gas (argon) occurs, we assume a uniform impact parameter distribution in the center-of-mass system (CMS) between the ion and the background gas. We first describe the simulations of DC sputtering, thereafter the simulations concerning high-power pulsed magnetron sputtering. The several interaction processes can be put into an abstract interaction model (pathway model, see [1]) that binds the interaction parameters together. A schematic drawing can be seen in Figure 3.

Figure 3: (a) Single (1 specie) pathway model, discussed in [1] and (b) Multiple (2 species) pathway model, discussed in [4].

4. Numerical Experiments

In the following, we discuss the numerical experiments, based on the influence of the particle flow via the BIAS voltage.

Here, the idea is to control the deposition rate with the BIAS voltage.

4.1. First Numerical Experiment: Delicate Geometry

In the first experiment, we study the influence of a delicate geometry to the DC-sputter process.

We apply an efficient particle-in-cell Monte Carlo simulation of an argon plasma in C++ (multiprocessor implementation via open MP). The simulations of arbitrary substrate geometries are performed, and ion sources and electrostatic boundary conditions are possible to implement.

We compute the spatial ion distribution between target and substrate (electrostatic biased) as well as spatial self-consistent electrostatic field configuration at equilibrium time (macroscopic time scale).

Next, we have the general setup in the reactor; see Figure 9.

For the computations, we studied the following objects as target; see Figure 10.

The parameters for the electric-field equations are given in the following Table 3.

Table 3: Parameters of the electrons.

In experiment A, we choose the and , the results are given in Figure 4, in experiment B, we choose and , the results are given in Figure 5, in experiment C, we choose and , and the results are given in Figure 6.

Figure 4: Target experiment with and .
Figure 5: Target experiment with and .
Figure 6: Target experiment with and .
4.1.1. A Phenomenological Study: Deposition Rate (Object A)

We measured the equilibrium deposition rate with respect to the BIAS voltage at the substrate. In order to estimate the deposition rate, we used the first 1500 time steps for coming to equilibrium and the following 1500 time steps for our measurement.

We apply two different particle depositions with drift of the species:  [m/s] and  [m/s].

We vary the BIAS between 1 [V] and 200 [V], but the variations did not affect the deposition rate of the particles, which are about (particles/timestep).

Remark 2. In the experiments, we show that the influence of the deposition rates at macroscopic time scales seems with respect to the underlying BIAS voltage.
In these experiments, we neglect the magnetic field in the Lorentz forces, and it seems that the deposition rates are independent of the BIAS voltage.

4.2. Second Numerical Experiment: Delicate Geometry (Concerning the Magnetic Field in the Lorentz Force on Each Particle)

In the second experiment, we have taken into account the full Lorentz force on each particle.

We simulate the paths of the sputtered atoms Ti and C and obtain at least a clue of the stoichiometric decomposition at the substrate.

We apply the full model (full Lorentz force) of an interaction between the sputtered particles and the background gas (Ar).

Due to relatively low-plasma densities, we consider our particle-in-cell (PIC) methods and study the plasma behavior by computing the trajectories of finite-size particles under the action of an external and self-consistent electric field defined in a grid of points. For the shown results, we use a computer cluster of Intel(R) Xeon(R) CPU X5472 at 3.00 GHz and a total memory of 64 GByte (due to large electrostatic and magnetostatic grids).

We solve the ion-electrostatic field feedback mechanism (self-regulating dynamic mechanism) and complex fields due to biased complex substrate geometries (mixed electrostatic boundary conditions) with the discussed methods.

By decoupling the ions and electrons (electrons behave like perfect fluid, due to lower mass), we could save computational time.

The Monte Carlo particle-in-cell simulations are only done for the transport of the ions, and we use the order of 10 mio. particles.

Based on the electric field, we have a transition from event-related MC to time-related MC (while computing the electrical influence to the ions).

Here, we study the influence of a delicate geometry with respect to the Lorentz forces of a DC-sputter process.

We have the general setup in the reactor; see Figure 11.

For the computations we studied the following objects as target: as seen in Figure 10.

The parameters for the electric-field equations are given in the Table 3.

In experiment A, we choose the Ar-Density =  mbar, .

The results are given in Figure 7.

Figure 7: Target experiment with Ar-Density = mbar, .

In experiment B, we choose the Ar-Density =  mbar, .

The results are given in Figure 8.

Figure 8: Target experiment with Ar-Density = mbar, .
Figure 9: The setup of the reactor for the experiments.
Figure 10: The following objects are used as targets for the experiments.
Figure 11: The general setup of the reactor for the experiments.

In experiment C, we choose a realistic test geometry by cooperation partner M. Balzer (FEM, Schwäbisch Gmünd, Germany) with the parameters: Ar-Density =  mbar, , and distance = 60 mm.

The results are given in Figure 12.

Figure 12: Target experiment (realistic test geometry) with Ar-Density = mbar, , and Distance = 60 mm.

In experiment D, we test the influence of a nonplanar substrates with the parameters: Ar-Density =  mbar and different distances.

The results are given in Figure 13.

Figure 13: Target experiment (realistic nonplanar substrates) with Ar-Density = mbar and different distances.

Remark 3. In the experiments, we show that the influence of the deposition rates at macroscopic time scales seems with respect to the underlying BIAS voltage. For the full model (with magnetic field in the Lorentz force), we obtain a dependency of the BIAS voltage to the deposition rate. In the arbitrary substrate geometries, we obtain the best stoichiometric composition at BIAS voltage of about −30 V and moderate target-substrate distances around 8 cm and a moderate Ar-Density of about  mbar. Further, we found out that an important quantity for nonplanar substrates is the ratio of the width and the depth of an inlet. So it makes sense to have at least small inlets to obtain a homogeneous deposition on the target.

5. Conclusions and Discussions

We present a novel model to consider DC and HIPIMS sputter processes in real gas regimes. We extend the underlying transport model by a Maxwell equation to model the electromagnetic field. The electrostatic field based on an additional BIAS voltage allows control of the ionized particles. We can improve the deposition rate in delicate geometrical zones. Mathematically, we solve coupled transport and Maxwell equations, while we apply the optimal solvers for each separate equations. Numerical results are presented and compared to real life applications, that allow to estimate an optimal BIAS voltage about −5 [V] to −30 [V].

5.1. Derivation of the Mean Relative Velocity

Within the framework of statistical mechanics, the mean value of an observable O can be computed via With describing canonical coordinates and canonical momenta of an -particle system, obey the Hamilton equation of motion. The probability distribution depends on the total Hamilton function of the system. Within the canonical ensemble, one has the following relation: With and the assumption of an ideal gas, the Hamilton function for the background gas is constructed only by the kinetic energies of the gas particles If , then the coordinate integration gives a volume factor in the numerator and denominator and therefore no contribution. The momentum integration can be done immediately and results in Gaussian integrals. The result for the mean relative velocity is therefore given by with , the reduced partition function (Maxwell distribution), and . By substituting and , one gets By using spherical coordinates with , , and , one gets The double integral on the right hand side can be evaluated, and its solution is given by After some simplification, the mean relative velocity reads whereby we made use of the scalar .

With , the final result for the mean relative velocity between projectiles probing into a monoatomic ideal gas is given by


  1. D. J. Christie, “Target material pathways model for high power pulsed magnetron sputtering,” Journal of Vacuum Science and Technology A, vol. 23, no. 2, pp. 330–335, 2005. View at Publisher · View at Google Scholar
  2. S. Longo, M. Capitelli, and K. Hassouni, “The coupling of PIC/MCC models of discharge plasmas with vibrational and electronic kinetics,” Journal De Physique, vol. 7, no. 4, pp. 271–281, 1997. View at Google Scholar
  3. D. J. Christie, “Particle-in-cell with Monte Carlo collisions gun code simulations of a surface-conversion H- ion source,” Communications in Computational Physics, vol. 4, no. 3, pp. 659–674, 2008. View at Google Scholar
  4. J. Geiser and R. Roehle, “Modeling and simulation for physical vapor deposition: multiscale model,” Journal of Convergence Information Technology, Korea, Article ID 1-075033, 2008. View at Google Scholar
  5. P. Eklund, Multifunctional Nanostructured Ti-Si-C Thin Films, vol. 1087 of Linköping studies in science and technology, Dissertation, Linköping University, 2007.
  6. S. Mahieu, G. Buyle, D. Depla, S. Heirwegh, P. Ghekiere, and R. De Gryse, “Monte Carlo simulation of the transport of atoms in DC magnetron sputtering,” Nuclear Instruments and Methods in Physics Research B, vol. 243, no. 2, pp. 313–319, 2006. View at Publisher · View at Google Scholar
  7. J. Geiser and S. Blankenburg, Monte Carlo Simulations Concerning Elastic Scattering with Application to DC and High Power Pulsed Magnetron Sputtering for Ti3 Si C2, Humboldt University of Berlin, Department of Mathematics, Berlin, Germany, 2009.