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

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.


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 SiH 4 /H 2 gas mixture during a PECVD process. Finally the conclusion summarizes the contents of the chapter.

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.  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 Si n H m oligomers (n ≤ 5). It presents a simulation of the growth of the films. Gorbachev et al. found that Si 2 H 5 and Si 3 H 7 strongly influence the growth of the film [11].
Valipa et al. [13] calculated the β reactivity of the SiH 3 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 CH 4 /H 2 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 SiH 4 /H 2 and CH 4 /H 2 gas mixtures during PECVD processes. The used numerical simulation techniques were MCS and MDS. To complete the studies, we used the fluid model [25].

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:

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.

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

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

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.

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): Let: Let x i and y i 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 r 1 and r 2 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 x i and y i are random numbers (x i ∈ {a, b} and y i ∈ {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 y i ≤ f(x i )). The number of favorable cases is N fav . 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 (x i , y i ) satisfying the condition: x i 2 + y i 2 ≤ 1. The function f(x) is equal to ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À x 2 p . 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).

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: a. 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 N pc particles (and/or elements). b. Choice of an initial configuration that responds to some physical and thermodynamic properties. The total or internal energy of the system is E i .
c. Infinitesimal random displacement of a particle (or element of the system) and calculation of the new internal energy of the system E f . 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 = E f -E i the fluctuation.
d. If ΔE ≤ 0; the new configuration is retained (favorable) and the different averages can be obtained; go to step (c).
e. If ΔE > 0; a random number ε is chosen such that 0 < ε < 1. Let the probability P r equal to: P r = exp. (ÀΔE/k B T) (where k B is the Boltzmann constant and T is the temperature).
f. If ε < P r , 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 P r . Figure 2 shows how to choose between the selected configurations. Let ε be a random number following a uniform law; If ε 1 ≤ P r the configuration is retained, and if ε 2 > P r the configuration is rejected.
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.

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 k B is called Boltzmann's constant. It is conventional to denote the quantity (k B T) À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: and 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.

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, surfaceinduced 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 r i belonging to the domain {0, 1} (or r i ∈ {0, 1} is available. This distribution law is uniform. To have a real random number x i belonging to the domain {a, b} (or x i ∈ a; b f g) according to a law of uniform distribution, we have: To have a real random number x i belonging to the domain {a, b} (or x i ∈ {a,b}) according to a formula (or law) of nonuniform distribution f(x), a histogram technique is used. Let N m be the number of intervals. If the mesh is regular (Figure 3): We define: We define the sequence: and the sequence: Hence each real random number r i belongs to the domain {0, 1} (where r i ∈ {0, 1}) (according to the uniform law); this number belongs to the domain {rx j-1 , rx j }. It corresponds to a random value x ran of the domain {x j-1 , x j }; 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 Δx i , or with tabular data f(x i ) 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.

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.
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.

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 SiH 4 /H 2 mixture used in PECVD process. Some paragraphs have been treated in previous works [21,24].

Description of the physical phenomenon
We use a MCS to study collisions and chemical reactions in gas phase of SiH 4 /H 2 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 (SiH 2 ), and silyl (SiH 3 ). The hydrogen consumption reactions SiH 4 + H ! SiH 3 + H 2 and SiH 3 + H ! SiH 2 + H 2 are found to play a central role in deciding the distribution of hydrogen [39].The plasma chemistry indicates that H atoms and SiH 3 radicals play an important role in the a-Si:H deposition process [40]. Experimentally, it is generally accepted that SiH 3 radicals dominate a-Si:H and μc-Si film growth from SiH 4 plasmas in the PECVD; it is the key precursor of a-Si:H deposition [41]. The proposed MCS allowed to get the ratio SiH 2 /SiH 3 and mean Schematic of a hybrid model of three modules used to study gas mixtures in the PECVD [33].
value of densities of species. It provides information on SiH 4 dissociation and on the production of SiH 3 , H, SiH 2 , and Si 2 H 6 and other important parameters.
The plasma in the PECVD reactor is weakly ionized. For our study, the mixture gas contains 22% of SiH 4 and 78% of H 2 ; 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. 10 8 cm À3 . The process is considered to be stationary. We take into account electrons and eight neutral species (SiH 4 , SiH 3 , SiH 2 , H, H 2 , Si 2 H 6 , Si 2 H 5 , SiH). Reactions taken into account include seven electron-neutral and 14 neutral-neutral reactions. Table 1 shows the 21 reactions and rate constants K reac . 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 T e .
Let K cons reac and K prod reac 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:

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 L x , L y , and L z ( 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.
Let n i be the density of neutral spice i (i = 1,…, 8). The first particle i is randomly chosen according to a probability of neutral species Pr sp,I (nonuniform discrete distribution) given by: The chosen particle takes randomly three components of space in cell r i (x i , y i , z i ) according to the normal distribution (nonuniform distribution). It takes also randomly three components of velocity v i (vx i , vy i , vz i ) according to Maxwell-Boltzmann distribution.

Treatment of elastic and inelastic collisions
Let n i and n j be the densities of species i and j in the gas and V ij 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 <s ij > 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 Pr col,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 (E a ) 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 K reac is the rate constant of the gas phase reaction.
The two colliding particles (e.g., the electron and SiH 4 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 Pr reac (i,j): where P K reac i; j ð Þ 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 K reac 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: a. 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.
b. Choices of random collision with a spice j.
c. 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.
d. At each step, we can note the different statistics.

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 N i of particles of each species is made. A first estimate of the sizes (L x , L y , L z ) 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 N i of particles and the sizes (L x , L y , L z ) of the elementary cell.
Let kp be the number of a species, kp = 1,…, 9. The minimal numbers Qnp(kp) and the sizes (L x , L y , L z ) 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 N col,m be the maximum number of elastic collisions per particle, and let N cycle be the number of cycles to average the simulation calculations.
For this MCS, the numerical chosen values are in Table 2.
For radicals (e.g., SiH 3 ), 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.

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 Si 2 H 6 , SiH, and Si 2 H 5 particles reaching the surface is negligible.
• Let N s,i and N s, H 2 be the densities of a species i and H 2 reaching the surface. The ratios N s,i /N s , H 2 are too small ( Table 3). Qnp(e) Qnp9 Table 2.
Used quantities and parameters in calculations for the gas temperature Tg = 520 K.
• Let N s,i be the density of a species i reaching the surface and N v,i the density of same species i in volume. The ratios N s,i /N v,i are too small ( Table 4); the surface effect is negligible.
• The reactions begin with the dissociation (consumption) of H 2 and SiH 4 by R5, R1, and R2 reactions.
• The production of SiH 3 is done by R8, and then there is production of SiH 2 by R12.
• The second important chemical reaction in the SiH 4 dissociation is R1: SiH 4 + e ! SiH 3 + H + e [24]. This result is compatible with that of Perkins et al. [51] and that of Doyle et al. [52].

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 SiH 4 /H 2 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  v, j 6.695 10 À6 7.965 10 À6 775 10 À6 Table 4. Ratios N s,i /N v,i of particles reaching the surface compared to volume.
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.

Author details
Fethi Khelfaoui* and Oumelkheir Babahani Physics Department, RPPS Laboratory, Faculty of Mathematics and Matter Sciences, Kasdi Merbah Ouargla University, Algeria *Address all correspondence to: fethi.khelfaoui@gmail.com; khelfaoui.fe@univ-ouargla.dz © 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.