Open access peer-reviewed chapter
By Fethi Khelfaoui and Oumelkheir Babahani
Submitted: March 25th 2019Reviewed: July 11th 2019Published: August 19th 2019
DOI: 10.5772/intechopen.88559
Downloaded: 503
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.
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.
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].
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:
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.
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.
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
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.
Calculation of the definite integral for a function f(x) of a single variable x on domain {a, b} has been proposed (Figure 1):
The integral of a function f(x).
Let:
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
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).
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:
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).
Choice of an initial configuration that responds to some physical and thermodynamic properties. The total or internal energy of the system is Ei.
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.
If ΔE ≤ 0; the new configuration is retained (favorable) and the different averages can be obtained; go to step (c).
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).
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.
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.
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:
or
and
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.
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.
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.
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].
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
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):
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.).
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.
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.
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.
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.
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].
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.
| Symbol | Reactions | Kreac (cm3/s) |
|---|---|---|
| R1 | SiH4 + e→SiH3 + H+e | k1 = 3 × 10−11 [42] |
| R2 | SiH4 + e→SiH2 + 2H + e | K2 = 1.5 × 10−10 [42] |
| R3 | SiH4 + e→SiH + H + H2 + e | K3 = 9.34 × 10−12 [42] |
| R4 | SiH4 + e→SiH2 + H2 + e | K4 = 7.19 × 10−12 [42] |
| R5 | H2 + e→2H + e | K5 = 4.49 × 10−12 [42] |
| R6 | Si2H6 + e→SiH3 + SiH2 + H + e | K6 = 3.72 × 10−10 [42] |
| R7 | Si2H6 + e→SiH4 + SiH2 +e | K7 = 1.1 × 1010× (1.(1./(1. + (0.63 × P)))) [43] |
| R8 | SiH4 + H→SiH3 + H2 | K8 = 2.8 × 10−11 × exp.(−1250/T) [44] |
| R9 | SiH4 + SiH2→Si2H6 | K9 = 1.1 × 1010 × (1.−(1./(1. + (0.63 × P)))) [43] |
| R10 | SiH3 + SiH3→SiH4 + SiH2 | K10 = 0.45 × 1.5 × 10−10 [44] |
| R11 | SiH4 + Si2H5→SiH3 + Si2H6 | K11 = 5 × 10−13 [42] |
| R12 | SiH3 + H→SiH2 + H2 | K12 = 2 × 10−11 [44] |
| R13 | SiH3 + Si2H6→SiH4 + Si2H5 | K13 = 4 × 10−10 × exp. (−2500/T) [44] |
| R14 | SiH2 + H→SiH + H2 | k14 = 2 × 10−11 [44] |
| R15 | Si2H6 + H→Si2H5 + H2 | K15 = 0.66 × 2.4 × 10−10 × exp. (−1250/T) [43] |
| R16 | Si2H6 + H→SiH4 + SiH3 | K16 = 0.34 × 2.4 × 10−10 × exp. (−1250/T) [44] |
| R17 | SiH + H2→SiH3 | K17 = 2 × 10−12 [43] |
| R18 | SiH2 + SiH3→Si2H5 | K18 = 3.77 × 10−13 [43] |
| R19 | SiH2 + H2→SiH4 | K19 = 3 × 10−12 × (1. + (1./1. + (0.03 × P))) [43] |
| R20 | 2SiH3→Si2H6 | K20 = 0.1 × 1.5 × 10−10 [43] |
| R21 | SiH4 + SiH→Si2H5 | K21 = (1.−(1./(1. + (0.33 × P)))) × (6.9 × 10−10) [43] |
List of gas phase reactions and corresponding rate constants [24].
Let
And chemical reaction for the production of A is as:
Rate production and consumption for any species A are taken as:
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.
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.
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
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
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:
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.
Choices of random collision with a spice j.
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.
At each step, we can note the different statistics.
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 collisions | Number of species Kp | Initial number of particles in cell | ||
|---|---|---|---|---|
| Lx (m) | 4.68 10−6 | 1 | Qnp(SiH4) | Qnp1 |
| Ly (m) | 4.68 10−6 | 2 | Qnp(SiH3) | 10 |
| Lz (m) | 20.0 10−3 | 3 | Qnp(SiH2) | 10 |
| 4 | Qnp(H) | 10 | ||
| Ncol,m | 500 | 5 | Qnp(H2) | Qnp5 |
| Ncycle internal cycle | 2000 | 6 | Qnp(Si2H6) | 10 |
| Ncycle external cycle | 200,000 | 7 | Qnp(SiH) | 10 |
| 8 | Qnp(Si2H5) | 10 | ||
| 9 | Qnp(e) | Qnp9 | ||
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).
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].
| Type | H2 | SiH4 | H | SiH3 | SiH2 |
|---|---|---|---|---|---|
| Ns,i/Ns, H2 | 1 | 0.23 | 1.67 10−4 | 8.60 10−5 | 9.86 10−6 |
Ratios Ns,i/Ns, H2 of particles reaching the surface compared to H2.
| Type | SiH4 | SiH3 | SiH2 |
|---|---|---|---|
| v, j | 6.695 10−6 | 7.965 10−6 | 775 10−6 |
Ratios Ns,i/Nv,i of particles reaching the surface compared to volume.
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.
© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
503total chapter downloads
Login to your personal dashboard for more detailed statistics on your publications.
Access personal reportingEdited by Pooneh Saidi
By Artur Chrobak, Grzegorz Ziółkowski and Dariusz Chrobak
Edited by Shaul Mordechai
We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.
More About Us