Open access peer-reviewed chapter

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase during Thin Film Deposition

Written By

Fethi Khelfaoui and Oumelkheir Babahani

Submitted: March 25th, 2019 Reviewed: July 11th, 2019 Published: August 19th, 2019

DOI: 10.5772/intechopen.88559

Chapter metrics overview

1,089 Chapter Downloads

View Full Metrics


Many physical phenomena can be modeled using Monte Carlo simulation (MCS) because it is a powerful tool to study thermodynamic properties. MCS can be used to simulate interactions between several particles or bodies in the presence of local or external fields. The main idea is to create a high number of different random configurations; statistics can be taken according to appropriate algorithms (Metropolis algorithm). In this chapter, we present basic techniques of MCS as the choice of potential, reaction rates, simulation cell, random configurations, and algorithms. We present some principal ideas of MCS used to study particle-particle collisions in the gas and in plasmas. Other MCS techniques are presented briefly. A numerical application is presented for collisions in gas phase during thin film deposition by plasma-enhanced chemical vapor deposition (PECVD) processes. Parameters and results of the simulation are studied according to a chosen reactor and mixture.


  • MCS
  • potential
  • reaction rates
  • collisions
  • thin film deposition

1. Introduction

In statistical physics only a few problems can be solved exactly. For complex problems, numerical methods can give exact results for problems that could only be solved in an approximate way. Numerical simulation can be a way to test the theory. The numerical results can be compared to the experimental results. The numerical simulation is placed between the fundamental and the experimental treatment; it has a quasi-experimental character (numerical experience). For problems of statistical physics, the most widely used simulation methods are the Monte Carlo method and the molecular dynamics method.

The first Monte Carlo simulation (MCS) was proposed by Metropolis et al. in 1953 [1]. The second Monte Carlo simulation was proposed by Wood and Parker in 1957 [2]. The obtained results were in good agreement with the experimental results of Bridgman [3] and those of Michels et al. [4]. In this method we attribute a series of initial positions chosen randomly to a system of N particles interacting through a defined potential. A sequence of particle configurations is generated by giving successive displacements to particles; we only retain configurations to ensure that the probability density is that of the chosen.

Molecular dynamics simulation (MDS) has been first introduced to simulate the behavior of fluids and solids at the molecular or atomic level. MDS was used for the first time by Alder and Wainwright in the late 1950s [5, 6] to study the interactions of hard spheres. The principle is the resolution of equations of motion for a hard sphere system in a simulation cell. The basic algorithm is Verlet’s algorithm [7].

In this chapter, we will present techniques of numerical simulations using the Monte Carlo method. We will present an application on the gas phase during plasma-enhanced chemical vapor deposition (PECVD) of thin films. The application concerns collisions between particles. Particles are in Brownian motion. Collisions, elastic or inelastic, are considered to be binary. Non-elastic collisions result in effective chemical reactions.

In Section 2, we cite some MCS and MDS works on PECVD processes. Section 3 presents general rules on numerical simulation methods. Section 4 presents how to simulate a physical problem using MCS? We present the Metropolis algorithm as a scheme to trait random configurations and different modules related to elaborate an MCS code. In Section 5, we apply the MCS on SiH4/H2 gas mixture during a PECVD process. Finally the conclusion summarizes the contents of the chapter.


2. Simulation works on the PECVD using MCS and MDS

The PECVD is the most widely used technique to produce hydrogenated amorphous silicon thin films (a-Si:H) for solar cells and for film transistors and electronic devices [8, 9]. Reactions during plasma deposition are complex and are not understood completely.

Gorbachev et al. [10, 11, 12] have developed a model that is based on chemical reactions and different processes in a PECVD reactor. The model takes into account the formation of SinHm oligomers (n ≤ 5). It presents a simulation of the growth of the films. Gorbachev et al. found that Si2H5 and Si3H7 strongly influence the growth of the film [11].

Valipa et al. [13] calculated the β reactivity of the SiH3 radical on a surface of a silicon lattice plane during the growth of a-Si:H using MDS. The mechanisms of physical and chemical interactions of low temperature plasmas with surfaces can be explored using MDS [14].

For a CH4/H2 mixture, Farouk et al. used the Monte Carlo method (PIC/MC); they calculated the ionization rate of the plasma and the deposition rate of the thin layer [15]. Rodgers et al. [16] have developed three-dimensional Monte Carlo simulations of diamond (100) surface CVD. Other works on MCS are in [17, 18, 19].

In our previous works [20, 21, 22, 23, 24], we were interested in the study of the gas phase and the interaction of plasmas with the surface, for SiH4/H2 and CH4/H2 gas mixtures during PECVD processes. The used numerical simulation techniques were MCS and MDS. To complete the studies, we used the fluid model [25].


3. General rules for numerical simulation methods

The starting point of numerical simulation is a physical phenomenon; its purpose is to obtain useful physical results. Between these two points, several steps can be identified. These steps are general and they are applicable for MCS. The steps can be summarized as follows:

3.1 Definition of the physical phenomenon and main hypothesis

The physical phenomenon must be defined by the description of the dominant domain of physics. The main assumptions and simplifying approximations are necessary to understand the physical phenomenon and the design of the first model.

3.2 Definition of the mathematical model

Mathematical model requires a mathematical formulation of the problem. It may be a problem of elements or discrete object or a problem of a continuous medium; it may be a spatiotemporal problem or frequency problem and may be a deterministic or probabilistic problem.

It would be interesting to know the mathematical equations that govern the phenomenon:

  • The forces between particles and elements

  • The potential interaction

  • The determination of a time scale

  • The determination of a length scale

  • Definition of constant magnitudes of motion and equilibrium magnitudes

  • Continuity equations, balance equations, transfer equations, etc.

3.3 Elaboration of simulation code

The MCS technique has been chosen for this work; knowing its basic algorithm is necessary for elaborating the simulation. This step requires some actions:

  • Validation of the model on simple cases

  • Simulation calculation on complex phenomena


4. Algorithms and techniques for MCS

The MCS is based on a probabilistic process with a random choice of configurations and samples of the situation of the physical system. The two pedagogical examples most cited in the literature are the integration of a single variable function and Ising’s model of spin. In the following subsection, we define the integration of a single variable function. We introduce the Ising model at the end of Section 4.2.2.

4.1 Integration of function of a single variable

Calculation of the definite integral for a function f(x) of a single variable x on domain {a, b} has been proposed (Figure 1):

Figure 1.

The integral of a function f(x).



Let xi and yi be real random numbers (i = 1, 2,…, N), and let H be a real number greater than the f(x) for x belonging to the domain {a, b} (or x ∈ {a, b}).

Let r1 and r2 be two random numbers belonging to the domain {0, 1} according to a uniform distribution law. Generators (e.g., Ran, RANDOM, RANDUM, or other IMSL mathematical libraries) of random numbers can be used:


where xi and yi are random numbers (xi ∈ {a, b} and yi ∈ {0, H}).

The Monte Carlo (MC) method is based on a probabilistic process. Let N be the total number of cases chosen (possible cases). It is necessary to count the number of favorable cases (or the number of points below the curve y = f(x)); let yi ≤ f(xi)). The number of favorable cases is Nfav. When N➔∞, the value I of the integral is [26]:


An example [26] is the calculation of the value π by calculating the integral I on a quarter circle of unit radius (R = 1.0). The pairs of random numbers (xi, yi) satisfying the condition: xi2 + yi2 ≤ 1. The function f(x) is equal to 1x2.

We take a = 0.0, b = 1.0, and H = 1.0.

For different values of N, we show that the numerical solution tends to π = 4I.

Although this integral is simple, it shows the strength and simplicity of the method. The technique can be generalized for the integration of multivariate functions.

We note that integration by the MC method is based on:

  • The choice of random configurations according to a uniform distribution law

  • Each configuration chosen is either favorable or unfavorable (the “or” is exclusive).

4.2 Principle of the MCS model

4.2.1 Calculation algorithm (Metropolis algorithm)

For statistical physics problems, the probabilistic choice of configurations is not always deterministic; the favorable and unfavorable cases are not exclusive. According to the Metropolis algorithm [26, 27], the steps of the simulation are:

  1. Choice of a simulation cell of adequate shape to the studied phenomena. The size of the simulation cell is related to a scale of length characteristic of the forces and interaction potential of the studied phenomenon. This cell may contain Npc particles (and/or elements).

  2. Choice of an initial configuration that responds to some physical and thermodynamic properties. The total or internal energy of the system is Ei.

  3. Infinitesimal random displacement of a particle (or element of the system) and calculation of the new internal energy of the system Ef. This displacement is related to the physical magnitudes: time scale and length scale. The physical system tends toward a minimization of the internal energy of the system with some fluctuation. Let ΔE = Ef-Ei the fluctuation.

  4. If ΔE ≤ 0; the new configuration is retained (favorable) and the different averages can be obtained; go to step (c).

  5. If ΔE > 0; a random number ε is chosen such that 0 < ε < 1. Let the probability Pr equal to: Pr = exp. (−ΔE/kBT) (where kB is the Boltzmann constant and T is the temperature).

  6. If ε < Pr, accept the move and in any case go back to step (c) for a new choice of an infinitesimal displacement (new configuration). Note that if such a trial move is rejected, the old configuration is again counted in the averaging with probability Pr.

Figure 2 shows how to choose between the selected configurations. Let ε be a random number following a uniform law; If ε1 ≤ Pr the configuration is retained, and if ε2 > Pr the configuration is rejected.

Figure 2.

Configuration choice according to Metropolis scheme.

Numerical simulation using the MC method is a very important tool for the study of static properties. The basic algorithm is based on probability notions. Understanding of the distribution function and/or interaction potentials is the heart of the calculation.

4.2.2 Thermodynamic quantities at equilibrium

In equilibrium statistical physics, the system has a certain probability that can be in any states. The probability of being in a state μ with energy H(μ) is given by the Boltzmann distribution P(μ):


where T is the absolute temperature and kB is called Boltzmann’s constant. It is conventional to denote the quantity (kBT)−1 by the symbol β. The normalizing factor Z, or partition function, is given by:


The average of a quantity Q fora system in equilibrium is:


The internal energy U, is given by:


which can be written in terms of a derivative of the partition function:


From thermodynamics we have expressions for the specific heat C, the entropy S, and the Helmholtz free energy F:








We can calculate other parameters affecting the system.

The Monte Carlo method is an excellent technique for estimating probabilities, and we can take advantage of this property in evaluating the results. The simplest and most popular model of a system of interacting variables in statistical physics is the Ising model. It consists of spins σi which are confined to the sites of a lattice and which may have only the values (+1) and (−1). These spins interact with their nearest neighbors on the lattice with interaction constant J; they can interact with an external magnetic field B coupling to the spins. The Hamiltonian H for this model is [26]:


The Ising model has been studied in one and two dimensions to obtain results of thermal properties, phase transition, and magnetic properties [26, 27, 28]. For chosen values of J and/or B, different steps may be taken for the calculations (simulation cell, initialization, configurations, boundary conditions, calculation algorithms). For any configuration, each spin takes the two possible directions. The detail of the calculation procedure is not the purpose of this chapter.

4.2.3 MCS module designs Simulation cell and initialization

We give a system of N particles (atoms, molecules, ions or particles) placed in a cell of fixed volume, generally of cubic form. The initial positions may, depending on the case, be distributed randomly according to a certain law (uniform or otherwise) or have a given symmetry. In a fluid, a gas, or a plasma, the particles may have random positions in general; in a solid or surface, with a crystal structure, the particles take ordered positions. The choice of random initial positions allows great freedom on the choice of the number of particles in the cell.

At the first step, the particles are given velocities that are generally selected to have a zero total momentum. If the system is in thermodynamic equilibrium, the initial velocities will be randomly chosen according to a Maxwell-Boltzmann law. In the general case, the velocity distribution is according to the problem dealt with. All other phase properties can be initialized to the particles; the main thing is the conservation of the total quantities of the system. Potentials of interaction

The particles interact with each other according to chosen interaction potentials. Since the interaction potentials are specific for each “numerical experiment,” the main part of the work consists in calculating the interaction energies for each proposed configuration.

The choice of interaction potentials is directly related to the mathematical formulation of the problem according to the state of the medium: fluid, gas, plasma, or solid. It can be Lennard-Jones potential, Coulomb potential, Debye potential, Morse potential, Stillinger-Weber potential, Born-Mayer potential, Moliere potential, or others. Boundary conditions

In general, two main boundary conditions are used: periodic boundary conditions (PBC) and minimum image convention (MIC) [29].

To minimize the surface effect, periodic boundary conditions (PBC) [30] are invariably imposed. The simulation cell is reproduced throughout the space to form an infinite mesh. We can simulate the properties of an infinite system. The particles that we follow are in the central cell; if a particle crosses a wall with a certain velocity, its image returns with the same velocity by the opposite wall. Under these conditions, the number of particles in the central cell, and consequently the density, is constant. These conditions also allow the conservation of the energy and the momentum of the system and do not introduce periodic effects (because of the interaction between particles).

According to the hypotheses and according to the geometry of the problem, other boundary conditions are proposed [26]. For example, in order to model thin films, the simulation cells are longitudinal and parallel to the film; one uses PBC in the directions parallel to the film. In the direction normal to the film, free edge boundary conditions can be used. In such cases, it may be appropriate to also include surface fields and surface interactions. In this way, one can study phenomena such as wetting, interface localization-delocalization transitions, surface-induced ordering and disordering, etc.

The core of the program includes calculating the potential energies of particle configuration and particle collisions. The interactions and collisions between particles can be elastic or inelastic; they can be binary or collective. For computation, the interaction energy of a particle with its neighbors is carried out by refocusing a base cell on the particle. This particle only interacts with particles in this region. This is called the “minimal image convention” (MIC) [1]. Sampling of random data

Generally, a RANDOM generator of real random numbers ri belonging to the domain {0, 1} (or ri ∈ {0, 1} is available. This distribution law is uniform.

To have a real random number xi belonging to the domain {a, b} (or xi ab) according to a law of uniform distribution, we have:


To have a real random number xi belonging to the domain {a, b} (or xi ∈ {a,b}) according to a formula (or law) of nonuniform distribution f(x), a histogram technique is used. Let Nm be the number of intervals. If the mesh is regular (Figure 3):

Figure 3.

Random number selection according to f (x) distribution.


We define:


We define the sequence:


and the sequence:


Hence each real random number ri belongs to the domain {0, 1} (where ri ∈ {0, 1}) (according to the uniform law); this number belongs to the domain {rxj-1, rxj}. It corresponds to a random value xran of the domain {xj-1, xj}; this number satisfies the formula (or the law) of nonuniform distribution f(x).

This technique can be generalized for a nonuniform distribution law f(x) with an irregular mesh Δxi, or with tabular data f(xi) with i = 1,…, m.

The technique can be generalized, too, for a discrete distribution law f(i) with i = 1,…, m.

In the literature, the reader can find simple algorithms for the choice of random numbers of some simple functions (Gaussian, etc.). Control of the evolution of the physical system

It is necessary to find some parameters allowing the control of the smooth course of the evolution of the system. We must look for the constants of movement. For example for an isolated system, we have the conservation of the total energy and the quantity of matter. Statistical calculations

By using the numerical simulation, it is possible to calculate many spatiotemporal quantities F(r,t). These quantities can be positions, speeds, kinetic moments, particle energies, concentrations, transport coefficients, etc. It would then be possible to calculate all other quantities related to F(r,t).

For the calculation of the averages, one can note the quantities on the space, on the time or on both. The histogram methods can be used. Static or dynamic distribution functions and spatial or temporal correlation functions can be calculated. It should be noted that the SMC is much more adequate for static properties because of the probabilistic choice of configurations.

Any calculated function or parameter F(r,t) can be used for another application in another calculation program.

4.3 Other large methods of Monte Carlo simulation

In the MCS model discussed extensively in this chapter, it’s more about collisions between particles. It’s particle-particle MCS or PP-MCS. In many problems of physics, the general idea is the same, but the applications and proposed models are numerous.

Other MCS models, named particle-in-cell MCS (PIC-MCS), are based on particle-cell interactions. In these last models, we also use a probabilistic choice of configurations and small variations in the state of the system (following the Metropolis algorithm); the interaction is between the particle with a cell, a mesh, or a drop. The parameters and variables of the cell, although local and instantaneous, are macroscopic. These parameters and variables can be thermodynamic, fluid, or electromagnetic. An example of the model based on PIC-MCS is described by Mattei et al. [31] for simulation of electromagnetic particle-in-cell collision in inductively coupled plasmas. Several works can be found in the literature on this same line of work. Other MCS models using particles may be considered. [32].

For statistical physics problem solving (such as thin film deposition problems), MCS models use experimental, numerical, or theoretical data from other methods and models. Models can be improved to hybrid models. In the hybrid models, connections between two modules can be realized. The first module is MCS; the second module is fluid, electromagnetic, or other. An example of a three-module hybrid model is presented by Mao and Bogaerts [33] to study gas mixtures in PECVD system. The three modules are MCS, fluid, and electromagnetic. The first module EM calculates the electromagnetic fields by solving Maxwell equations. These fields are used as inputs in the module MCS, where the electron density, electron temperature, electron energy distribution function, and electron impact reaction rates can be computed with a Monte Carlo procedure. Subsequently, the module fluid calculates densities and fluxes of the various plasma species (i.e., heavy particles and electrons) with continuity equations and the electrostatic field with Poisson’s equation. This electrostatic field is used as input again in the EM. This cycle is iterated until convergence. The schematic of the hybrid model is given in Figure 4.

Figure 4.

Schematic of a hybrid model of three modules used to study gas mixtures in the PECVD [33].

To solve statistical physics problems with evolutions as a function of time, kinetic models of MCS (kMCS) are used. Using kMCS, Battaile and Srolovitz [17] described kinetic phenomena of the diffusive motion of a single interstitial atom in a close-packed metal crystal. The motion of the interstitial atom is usually limited to two types: vibration of the atom around the center of the interstitial hole in which it resides and hops to nearest-neighbor interstitial sites. The atom can hop into any of the nearest-neighbor interstitial sites; it executes a random walk. In an MC simulation of this diffusion process, the new position of the interstitial atom is chosen at random from a list of the adjacent interstitial sites.

Other CVD and PECVD works on MCS are presented in Ref.s [15, 34, 35, 36, 37, 38]. They show how MCS methods can study properties of gas mixtures and properties of the growth of thin films.


5. Example of application: Monte Carlo simulation of a gas mixture in the PECVD

In this section, we present an example of PP-MCS of collisions and reactions in gas phase of SiH4/H2 mixture used in PECVD process. Some paragraphs have been treated in previous works [21, 24].

5.1 Description of the physical phenomenon

We use a MCS to study collisions and chemical reactions in gas phase of SiH4/H2 mixture used in the PECVD process. In this phase, important reactions have been identified that contribute to the production and the consumption of hydrogen (H), silylene (SiH2), and silyl (SiH3). The hydrogen consumption reactions SiH4 + H → SiH3 + H2 and SiH3 + H → SiH2 + H2 are found to play a central role in deciding the distribution of hydrogen [39].The plasma chemistry indicates that H atoms and SiH3 radicals play an important role in the a-Si:H deposition process [40]. Experimentally, it is generally accepted that SiH3 radicals dominate a-Si:H and μc-Si film growth from SiH4 plasmas in the PECVD; it is the key precursor of a-Si:H deposition [41]. The proposed MCS allowed to get the ratio SiH2/SiH3 and mean value of densities of species. It provides information on SiH4 dissociation and on the production of SiH3, H, SiH2, and Si2H6 and other important parameters.

The plasma in the PECVD reactor is weakly ionized. For our study, the mixture gas contains 22% of SiH4 and 78% of H2; the pressure is 100 mtorr, the temperature of the gas ranges from 373 to 723 K, the electron temperature is about 2.5 eV, and the electron density is 3. 108 cm−3. The process is considered to be stationary. We take into account electrons and eight neutral species (SiH4, SiH3, SiH2, H, H2, Si2H6, Si2H5, SiH). Reactions taken into account include seven electron-neutral and 14 neutral-neutral reactions. Table 1 shows the 21 reactions and rate constants Kreac. At low temperature, the neutrals interact occasionally with each other and move under the effect of thermal agitation; their velocity distribution function is Maxwell-Boltzmann distribution. Electrons have the mean velocity with kinetic energy Te.

SymbolReactionsKreac (cm3/s)
R1SiH4 + e→SiH3 + H+ek1 = 3 × 10−11 [42]
R2SiH4 + e→SiH2 + 2H + eK2 = 1.5 × 10−10 [42]
R3SiH4 + e→SiH + H + H2 + eK3 = 9.34 × 10−12 [42]
R4SiH4 + e→SiH2 + H2 + eK4 = 7.19 × 10−12 [42]
R5H2 + e→2H + eK5 = 4.49 × 10−12 [42]
R6Si2H6 + e→SiH3 + SiH2 + H + eK6 = 3.72 × 10−10 [42]
R7Si2H6 + e→SiH4 + SiH2 +eK7 = 1.1 × 1010× (1.(1./(1. + (0.63 × P)))) [43]
R8SiH4 + H→SiH3 + H2K8 = 2.8 × 10−11 × exp.(−1250/T) [44]
R9SiH4 + SiH2→Si2H6K9 = 1.1 × 1010 × (1.−(1./(1. + (0.63 × P)))) [43]
R10SiH3 + SiH3→SiH4 + SiH2K10 = 0.45 × 1.5 × 10−10 [44]
R11SiH4 + Si2H5→SiH3 + Si2H6K11 = 5 × 10−13 [42]
R12SiH3 + H→SiH2 + H2K12 = 2 × 10−11 [44]
R13SiH3 + Si2H6→SiH4 + Si2H5K13 = 4 × 10−10 × exp. (−2500/T) [44]
R14SiH2 + H→SiH + H2k14 = 2 × 10−11 [44]
R15Si2H6 + H→Si2H5 + H2K15 = 0.66 × 2.4 × 10−10 × exp. (−1250/T) [43]
R16Si2H6 + H→SiH4 + SiH3K16 = 0.34 × 2.4 × 10−10 × exp. (−1250/T) [44]
R17SiH + H2→SiH3K17 = 2 × 10−12 [43]
R18SiH2 + SiH3→Si2H5K18 = 3.77 × 10−13 [43]
R19SiH2 + H2→SiH4K19 = 3 × 10−12 × (1. + (1./1. + (0.03 × P))) [43]
R202SiH3→Si2H6K20 = 0.1 × 1.5 × 10−10 [43]
R21SiH4 + SiH→Si2H5K21 = (1.−(1./(1. + (0.33 × P)))) × (6.9 × 10−10) [43]

Table 1.

List of gas phase reactions and corresponding rate constants [24].

Let Kreaccons and Kreacprod be the rate constants of the consumption and the production of species A. The chemical reaction for the consumption of A is as:


And chemical reaction for the production of A is as:


Rate production and consumption for any species A are taken as:


5.2 Description of Monte Carlo simulation technique

5.2.1 Simulation cell and phase coordinates

The MCS is based on binary collisions at the microscopic level. Elastic collisions are between all particles, and inelastic collisions (or effective collisions) are those that result in a chemical reaction. A chemical reaction needs a collision involving at least two particles (atoms, ions, electrons, or molecules). According to kinetic theory, gases consist of particles in random motion. These particles are uniformly distributed in a cell which has a parallelepiped form of sizes Lx, Ly, and Lz (Figure 5). These particles move in a straight line until they collide with other particles or the walls of their container. Dimensions and volume of Monte Carlo cell must take into consideration the mean free path of species.

Figure 5.

Form of the simulation cell.

Let ni be the density of neutral spice i (i = 1,…, 8). The first particle i is randomly chosen according to a probability of neutral species Prsp,I (nonuniform discrete distribution) given by:


The chosen particle takes randomly three components of space in cell ri(xi, yi, zi) according to the normal distribution (nonuniform distribution). It takes also randomly three components of velocity vi (vxi, vyi, vzi) according to Maxwell-Boltzmann distribution.

5.2.2 Treatment of elastic and inelastic collisions

Let ni and nj be the densities of species i and j in the gas and Vij the relative velocity between the two species i and j.

According to the kinetic theory of gases, we have for an incident particle i on a target particle j the average collision frequency νij as:


where <sij> is the cross section of the particle j.

The mean free path ι> of species i is:


The time between two collisions τij is then:


For chemical effective reactions (inelastic collisions) between two reactive species i and j giving products i’ and j’, the rate constant reaction verifies [45]:


General rules of collision theory are applied:

  • The new velocities of the colliding particles are calculated using conservation of energy and momentum for elastic collisions.

  • Conservation of total energy as isolated system.

  • Movement of the center of mass and relative motion around the center of mass.

The reader can refer to some fundamental physics books that deal with general notions of collisions and corresponding parameters [45, 46, 47, 48].

The plasma in the PECVD reactor is weakly ionized. At low temperature, particles interact occasionally with each other and move under the effect of thermal agitation. In reality, only a small fraction of collisions are effective (result in a chemical reaction) [21].

In our MCS, after traveling a random walk given by a Gaussian distribution, the first chosen particle collides with a second particle (molecule, atom, radical, or electron). The last particle j is randomly chosen according to a (i-j) collision probability Prcol,j (nonuniform discrete distribution) given by:


where νij is the neutral-neutral or electron-neutral collisional frequency. The collision theory indicates that the collision between molecules can provide the energy needed to break the necessary bonds so that new bonds can be formed [49]. Particles must have sufficient energy to initiate the reaction (activation energy), so the two chosen particles must have kinetic energy equal to or greater than the barrier energy (Ea) of a gas phase reaction. The difference between the kinetic energy of the two particles and the activation energy define the kind of collision (effective or not effective).

The activation energy is given by:


where the pre-exponential factor is assumed to be the collision frequency factor and Kreac is the rate constant of the gas phase reaction.

The two colliding particles (e.g., the electron and SiH4 molecule) can interact by several reactions (R1, R2, R3, and R4 in Table 1); we choose randomly one of gas phase reactions occurring according to a, nonuniform discrete distribution reaction probability Prreac (i,j):


where Kreacij is the sum of all rate constants of possible reactions between i and j.

All chemical systems go naturally toward states of minimum Gibbs free energy [21, 24]. A chemical reaction tends to occur in the direction of lower Gibbs free energy. To determine the direction of the reaction that is taking place, we use the old and new values of Kreac and the equilibrium constant with reactants and product concentrations. Each set of binary collisions can be related or converted into time. As cited in section (a), Table 1 gives gas phase reactions and corresponding rate constants used in this MCS.

To continue the simulation, after the elastic collision, particle i takes new values of components velocity and new mean free path; mean free path is taken from a normal (nonuniform) distribution (Gaussian distribution). If the collision is inelastic, we have to take a new particle.

From Metropolis algorithm, the scheme of this MCS is as follows:

  1. Choices of particle of spice i with random position, velocity, and mean free path; periodic boundary conditions are used to keep particles in the elementary cell.

  2. Choices of random collision with a spice j.

  3. Study of collision type (elastic, inelastic). If the collision is elastic the particle i move with a new velocity and mean free path, and we return to step (b). If the collision is inelastic particles i and j give new particles i’ and j’, according to Metropolis scheme, and we return to step (a) or (b). Periodic boundary conditions are used to keep particles in the elementary cell.

  4. At each step, we can note the different statistics.

5.2.3 The choice of simulation parameters

Once the species are selected for the simulation model, an estimate of species densities should be made. Following the model of interaction and collisions between particles (binary, collective, etc.), a first choice of the minimum number Ni of particles of each species is made. A first estimate of the sizes (Lx, Ly, Lz) of the elementary cell is made.

The study of the types of interaction potentials and the calculation of the approximate values of the force ranges, the kinetic energies, the internal energies, and the energies of activation make it possible to correct the minimal numbers Ni of particles and the sizes (Lx, Ly, Lz) of the elementary cell.

Let kp be the number of a species, kp = 1,…, 9. The minimal numbers Qnp(kp) and the sizes (Lx, Ly, Lz) have to be discussed for statistical calculations.

For numerical programming, according to the programming language used and according to the size (or the computational capacity) of the computer, it is necessary to find a judicious choice of the tables of integer or real values and which values would be useful to save all during simulation. Let Ncol,m be the maximum number of elastic collisions per particle, and let Ncycle be the number of cycles to average the simulation calculations.

For this MCS, the numerical chosen values are in Table 2.

Cell dimensions and steps for collisionsNumber of species KpInitial number of particles in cell
Lx (m)4.68 10−61Qnp(SiH4)Qnp1
Ly (m)4.68 10−62Qnp(SiH3)10
Lz (m)20.0 10−33Qnp(SiH2)10
Ncycle internal cycle20006Qnp(Si2H6)10
Ncycle external cycle200,0007Qnp(SiH)10

Table 2.

Used quantities and parameters in calculations for the gas temperature Tg = 520 K.

For radicals (e.g., SiH3), particle numbers Qnp(k) are very small; we take Qnp(k) = 10. These numbers cannot take value 1 or 0, even if a species k is in trace form in the gas. The value 0 for a species k means that any other species k’ does not make a collision with the species k; and the value 1 means that we have no collisions between particles of the same species in the cell.

Qnp1, Qnp5, and Qnp9 are calculated from the volume of cell, the pressure, the temperature, and the total number of particles in the cell (Qnp1 = 0.81187824 * 109; Qnp5 = 0.20296956 * 109; Qnp9 = 131).

5.2.4 Calculation of statistical properties and some results of the calculations

As we have chosen a stationary regime, we must reach the values and properties at equilibrium. The results of the simulation show this trend. In MCS, averaged values, distribution functions, autocorrelation functions, and correlation functions can be calculated. To ensure rapid convergence of calculations, it would be useful to look for statistically symmetric (or stationary or unsteady) parameters [26, 50].

As an example for our MCS calculation, we have:

  • The number of Si2H6, SiH, and Si2H5 particles reaching the surface is negligible.

  • Let Ns,i and Ns, H2 be the densities of a species i and H2 reaching the surface. The ratios Ns,i/Ns, H2 are too small (Table 3).

  • Let Ns,i be the density of a species i reaching the surface and Nv,i the density of same species i in volume. The ratios Ns,i/Nv,i are too small (Table 4); the surface effect is negligible.

  • The reactions begin with the dissociation (consumption) of H2 and SiH4 by R5, R1, and R2 reactions.

  • The production of SiH3 is done by R8, and then there is production of SiH2 by R12.

  • The reaction R2: SiH4 + e → SiH2 + 2H + e plays the central role in SiH4 dissociation by electron impact [24]. This result is compatible with [39].

  • The second important chemical reaction in the SiH4 dissociation is R1: SiH4 + e → SiH3 + H + e [24]. This result is compatible with that of Perkins et al. [51] and that of Doyle et al. [52].

Ns,i/Ns, H210.231.67 10−48.60 10−59.86 10−6

Table 3.

Ratios Ns,i/Ns, H2 of particles reaching the surface compared to H2.

v, j6.695 10−67.965 10−6775 10−6

Table 4.

Ratios Ns,i/Nv,i of particles reaching the surface compared to volume.


6. Conclusions

MCS is a widely used method in statistical physics to study thermodynamic, structural, or phase properties. It is based on random and probabilistic processes. The purpose of this chapter is to present the technique for general use in physics for the study of thin film deposition problems. The technique can be generalized to other fields of science: biology, economics, transportation, and social sciences.

We started by presenting general rules for numerical simulation methods. Metropolis algorithm has been considered as the basic algorithm. After, we presented the different steps for the realization of a MCS code. We chose the particle-particle model MCS (PP-MCS) to explain the different steps and procedures to be applied in the deposition of thin layers by PECVD processes. We have shown that this technique can be generalized to the particle-in-cell MCS (PIC-MCS) case or kinetic MCS (kMCS), as it can be joined with other modules to give hybrid models. It is important to know how to choose random configurations from the laws or probability distributions in the system.

A numerical application is presented for collisions in a SiH4/H2 gas mixture in the PECVD process. A preliminary work of determination of the chemical reactions between molecules and radicals is made. A choice of the simulation cell is made, and the definition of the probabilities of the collisions between peers is made. The Metropolis algorithm makes it possible to follow the various elastic and inelastic collisions; it also makes it possible to make the statistics of the interactions with the surface. The results are compatible with [39, 51, 52].

Other questions may be asked to account for molecular ions, surface and volume kinetics, or thin film formation. The techniques and different models of the MCS (PP-MCS, MCS-PIC, kMCS) allow taking care of these questions.

The interconnection of the MCS with other models (MDS, hybrid model, fluid model, electromagnetic model, etc.) would allow answering more questions. The methods can be applied to other specialties than the physical sciences.


  1. 1. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. The Journal of Chemical Physics. 1953;21:1087-1092. DOI: 10.1063/1.1699114
  2. 2. Wood WW, Parker FR. Monte Carlo equation of state of molecules interacting with the Lennard-Jones potential. The Journal of Chemical Physics. 1957;27:720-733. DOI: 10.1063/1.1743822
  3. 3. Bridgman PW. Compressibilities and electrical resistance under pressure, with special reference to intermetallic compounds. Proceedings of the American Academy of Arts and Sciences. 1935;70:285-317. DOI: 10.2307/20023138
  4. 4. Michels A, Wijker H, Wijker HK. Isotherms of argon between 0°C and 150°C and pressures up to 2900 atmospheres. Physica. 1949;15:627-632. DOI: 10.1016/0031-8914(49)90119-6
  5. 5. Alder BJ, Wainwright TE. Phase transition for a hard sphere system. The Journal of Chemical Physics. 1957;27:1208. DOI: 10.1063/1.1743957
  6. 6. Alder BJ, Wainwright TE. Studies in molecular dynamics. I. General method. The Journal of Chemical Physics. 1959;31:459-466. DOI: 10.1063/1.1730376
  7. 7. Verlet L. Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physics Review. 1967;159:98-103. DOI: 10.1103/PhysRev.159.98
  8. 8. Wong J, Kherani NP, Zukotynski S. Monte Carlo modeling of the dc saddle field plasma: Discharge characteristics of N2 and SiH4. Journal of Applied Physics. 2007;101:013308-013314. DOI: 10.1063/1.2409566
  9. 9. Matsuda A, Nomoto K, Takeuchi Y, Suzuki A, Yuuki A, Perrin J. Temperature dependence of the sticking and loss probabilities of silyl radicals on hydrogenated amorphous silicon. Surface Science. 1990;227:50-56. DOI: 10.1016/0039-6028(90)90390-T
  10. 10. Gorbachev YE, Zetevakhin MA, Krzhizhanovskaya VV, Shveïgert VA. Special features of the growth of hydrogenated amorphous silicon in PECVD reactors. Journal of Technical Physics. 2000;45:1032-1041. DOI: 10.1134/1.1307013
  11. 11. Gorbachev YE. Effect of oligomers on the growth of amorphous silicon films in a PECVD reactor. Technical Physics. 2006;51:733-739. DOI: 10.1134/S1063784206060089
  12. 12. Gorbachev YE, Zetevakhin MA, Kaganovich ID. Simulation of the growth of hydrogenated amorphous silicon films from an RF discharge plasma. Journal of Technical Physics. 1996;41:1247-1258
  13. 13. Valipa MS, Aydil ES, Maroudas D. Atomistic calculation of the SiH3 surface reactivity during plasma deposition of amorphous silicon thin films. Surface Science. 2004;572:L339-L347. DOI: 10.1016/j.susc.2004.08.029
  14. 14. Graves DB, Brault P. Molecular dynamics for low temperature plasma-surface interaction studies. Journal of Physics D. 2009;42:1-27. DOI: 10.1088/0022-3727/42/19/194011
  15. 15. Pandey SC, Singh T, Maroudas D. Kinetic Monte Carlo simulations of surface growth during plasma deposition of silicon thin films. The Journal of Chemical Physics. 2009;131:1-12. DOI: 10.1063/1.3152846
  16. 16. Rodgers WJ, May PW, Allan NL, Harvey JN. Three-dimensional kinetic Monte Carlo simulations of diamond chemical vapor deposition. The Journal of Chemical Physics. 2015;142:214707. DOI: 10.1063/1.4921540
  17. 17. Battaile CC, Srolovitz DJ. Kinetic Monte Carlo simulation of chemical vapor deposition. Annual Review of Materials Research. 2002;32:297-319. DOI: 10.1146/annurev.matsci.32.012102.110247
  18. 18. YangY G, ZhouX W, Johnson RA, Wadlty HNG. Monte Carlo simulation of hyperthermal physical vapor deposition. Acta Materialia. 2001;49:3321-3332. DOI: 10.1016/S1359-6454(01)00139-2
  19. 19. Chen Z-Y, Zhu Y, Chen S-H, Qiu Z-R, Jiang S-J. The kinetic process of non-smooth substrate thin film growth via parallel Monte Carlo method. Applied Surface Science. 2011;257:6102-6106. DOI: 10.1016/j.apsusc.2011.02.004
  20. 20. Babahani O, Khelfaoui F, Meftah MT. Analytical calculation of site and surface reaction probabilities of SiHx radicals in PECVD process. The European Physical Journal Applied Physics. 2013;62:10301-10307. DOI: 10.1051/epjap/2013120345
  21. 21. Babahani O. Simulation numérique par la méthode de Monte Carlo de la déposition de couches minces par procédés CVD [thesis]. Ouargla: Kasdi Merbah Ouargla University; 2013
  22. 22. Kebaili HO, Babahani O, Khelfaoui F. Simulation par la dynamique moléculaire de l'interaction plasma-surface lors de la croissance de couches minces a-Si:H par procédés PECVD. Annales des Sciences et Technologie. 2014;6:165-171
  23. 23. Babahani O, Khelfaoui F. Simulation Monte Carlo de réactions chimiques dans le volume d'un réacteur PECVD lors de déposition d'une couche mince a-C:H. Annales des Sciences et Technologie. 2014;6:156-164
  24. 24. Babahani O, Hadjadj S, Khelfaoui F, Kebaili HO, Lemkeddem S. Monte Carlo simulation of chemical reactions in plasma enhanced chemical vapor deposition: From microscopic view to macroscopic results. Silicon. 2019;11:1267-1274. DOI: 10.1007/s12633-018-9916-y
  25. 25. Babahani O, Khelfaoui F. Calcul des concentrations de molécules et de radicaux lors de déposition de couches minces a-Si:H par procédés PECVD. Annales des Sciences et Technologie. 2012;4:115-120
  26. 26. Landau DP, Binder K. A Guide to Monte Carlo Simulations in Statistical Physics. Cambridge University Press; 2009. DOI: 10.1017/CBO9780511994944
  27. 27. Murthy KPN. An Introduction to Monte Carlo Simulations in Statistical Physics. Tamilnadu, India: Indira Gandhi Centre for Atomic Research; 2003
  28. 28. Gould H, Tobochnik J. An Introduction to Computer Simulation Methods: Applications to Physical Systems, Part 2. Addisson-Wesley Publishing Company; 2006. ISBN-10: 0805377581
  29. 29. Hansen JP, Mc Donald IR. Theory of Simple Liquids. London: Academic Press; 1976. 428 p. ISBN: 9780080455075
  30. 30. Born M, Von Karman T. Reappraising 1907 Einstein’s model of specific heat. Physikalishce Zeitschrift. 1912;13:297-309
  31. 31. Mattei S, Nishida K, Onai M, Lettry J, Tran MQ, Hatayama A. A fully-implicit particle-In-cell Monte Carlo collision code for the simulation of inductively coupled plasmas. Journal of Computational Physics. 2017;350:891-906. DOI: 10.1016/
  32. 32. Hockney RW, Eastwood JW. Computer Simulation Using Particles. IOP Publishing Ltd; 1988
  33. 33. Mao M, Bogaerts A. Investigating the plasma chemistry for the synthesis of carbon nanotubes/nanofibers in an inductively coupled plasma enhanced CVD system: Effect of different gas mixtures. Journal of Physics D: Applied Physics. 2010;43:205201. DOI: 10.1088/0022-3727/43/20/205201
  34. 34. Battaile CC, Srolovitz DJ. A kinetic Monte Carlo method for the atomic-scale simulation of chemical vapor deposition: Application to diamond. Journal of Applied Physics. 1997;82:6293. DOI: 10.1063/1.366532
  35. 35. Lothar K, Kuhn Frank M, Olaf D. Kinetic Monte Carlo simulations of surface reactions on supported nanoparticles: A novel approach and computer code. The Journal of Chemical Physics. 2015;143:044108. DOI: 10.1063/1.4926924
  36. 36. May PW, Harvey JN, Allan NL, Richley JC, Mankelevich Yu A. Simulations of chemical vapor deposition diamond film growth using a kinetic Monte Carlo model and two dimensional models of microwave plasma and hot filament chemical vapor deposition reactors. Journal of Applied Physics. 2010;108:114909. DOI: 10.1063/1.3516498
  37. 37. Harris SJ. Mechanism for diamond growth from methyl radicals. Applied Physics Letters. 1990;56:2298. DOI: 10.1063/1.102946
  38. 38. May PW, Harvey JN, Allan NL, Richley JC, Mankelevich Yu A. Simulations of chemical vapor deposition diamond film growth using a kinetic Monte Carlo model. Journal of Applied Physics. 2010;108:014905. DOI: 10.1063/1.3437647
  39. 39. Aman-ur-Rehman, Kwon HC, Park WT, Lee JK. A study of the role of various reactions on the density distribution of hydrogen, silylene, and silyl in SiH4/H2 plasma discharges. Physics of Plasmas. 2011;18:093502. DOI: 10.1063/1.3630933
  40. 40. Moravej M, Babayan SE, Nowling GR, Yang X, Hicks RF. Plasma enhanced chemical vapour deposition of hydrogenated amorphous silicon at atmospheric pressure. Science and Technology. 2004;13:8-14
  41. 41. Hoefnagels JPM, Barrell Y, Kessels WMM, Van de Sanden MCM. Inflow and shock formation in supersonic, rarefied plasma expansions. Journal of Applied Physics. 2004;96:4094
  42. 42. Gorbachev YE, Zatevakhin MA, Krzhizhanovskaya VV, Shveigert VA. Special features of the growth of hydrogenated amorphous silicon in PECVD reactors. Technical Physics. 2000;45:1032-1041. DOI: 10.1134/1.1307013
  43. 43. Perrin J, Leroy O, Bordage MC. Cross-sections, rate constants and transport coefficients in silane plasma chemistry. Contributions to Plasma Physics. 1996;36:3-49. DOI: 10.1002/ctpp.2150360102
  44. 44. Amanatides E, Stamou S, Mataras D. Gas phase and surface kinetics in plasma enhanced chemical vapor deposition of microcrystalline silicon: The combined effect of RF power and hydrogen dilution. Journal of Applied Physics. 2001;90:5786-5798. DOI: 10.1063/1.1413241
  45. 45. Moison M, Pelletier J. Physique Des Plasmas Collisionnels: Application Aux décharges Haute fréquence. France: EDP Sciences; 2006
  46. 46. Delcroix J-L, Bers A. Physique Des Plasmas 1. Paris: Inter Editions/CNRS Editions; 1994
  47. 47. Atkins PW, de Paula J. Chimie Physique. 4th ed. De Boeck; 2015. 973 p. ISBN-10: 2804166511
  48. 48. Atkins PW. Eléments de Chimie Physique. De Boeck Université; 1998. 512 p. ISBN-10: 2744500100
  49. 49. Moore JT, Hren C, Mikulecky PJ. U Can: Chemistry I for Dummies. Wiley; 2015. 456 p. ISBN: 978-1-119-07940-8
  50. 50. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Oxford: Clarendon Press; 1987
  51. 51. Perkins GGA, Austin ER, Lampe FW. The 147-nm photolysis of monosilane. Journal of the American Chemical Society. 1979;101:1109-1115
  52. 52. Doyle JR, Doughty DA, Gallagher A. Silane dissociation products in deposition discharges. Journal of Applied Physics. 1990;68:4375

Written By

Fethi Khelfaoui and Oumelkheir Babahani

Submitted: March 25th, 2019 Reviewed: July 11th, 2019 Published: August 19th, 2019