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

Abstract

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.

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 SiH₄/H₂ 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 Si_nH_m oligomers (n ≤ 5). It presents a simulation of the growth of the films. Gorbachev et al. found that Si₂H₅ and Si₃H₇ strongly influence the growth of the film [11].

Valipa et al. [13] calculated the β reactivity of the SiH₃ 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₄/H₂ 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₄/H₂ and CH₄/H₂ 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:

$$I=\underset{a}{\overset{b}{\int }}f\left(x\right)\mathit{dx}$$E1

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₁ and r₂ 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:

$$x_i=a+r_1\left(b-a\right)\text{and}{y}_i=0+r_2\left(H-0\right)$$E2

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]:

$$I=\frac{N_{\mathit{fav}}}{N}\left(b-a\right)\left(H-0.0\right)=\frac{N_{\mathit{fav}}}{N}\left(b-a\right)H$$E3

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² + y_i² ≤ 1. The function f(x) is equal to $\sqrt{1-x^2}$.

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

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

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 P_r equal to: P_r = exp. (−ΔE/k_BT) (where k_B is the Boltzmann constant and T is the temperature).

6. 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 ε₁ ≤ P_r the configuration is retained, and if ε₂ > P_r 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(μ):

$$P_{\mu }=\frac{exp\left(-H\left(\mu \right)/k_{B}T\right)}{Z}$$E4

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

$$Z=\sum _{\mu }exp\left(-H\left(\mu \right)/k_{B}T\right)=\sum _{\mu }exp\left(-\mathit{\beta H}\left(\mu \right)\right)$$E5

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

$$<Q>=\sum _{\mu }{Q}_{\mu }{P}_{\mu }=\frac{1}{Z}\sum _{\mu }{Q}_{\mu }exp\left(-\mathit{\beta H}\left(\mu \right)\right)$$E6

The internal energy U, is given by:

$$U=\frac{1}{Z}\sum _{\mu }H\left(\mu \right)exp\left(-\mathit{\beta H}\left(\mu \right)\right)$$E7

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

$$U=\frac{1}{Z}\frac{\partial Z}{\partial \beta }=-\frac{\partial logZ}{\partial \beta }$$E8

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

$$C=\frac{\partial U}{\partial T}=-k_B{\beta }^2\frac{\partial U}{\partial \beta }=-k_B{\beta }^2\frac{{\partial }^{2}logZ}{\partial {\beta }^2}$$E9

or

$$C=T\frac{\partial S}{\partial T}=-\beta \frac{\partial S}{\partial \beta }$$E10

and

$$S=-k_B\beta \frac{\partial logZ}{\partial \beta }+k_{B}logZ$$E11

and

$$F=U-\mathit{TS}=-k_{B}logZ$$E12

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]:

$$H=-J\sum _{i,j}{\sigma }_i{\sigma }_j-B\sum _i{\sigma }_i$$E13

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

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

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

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

4.2.3.4 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 $\in \left\{a,b\right\}$) according to a law of uniform distribution, we have:

$$x_i=a+r_i\left(b-a\right)$$E14

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):

Figure 3.

Random number selection according to f (x) distribution.

$$\mathrm{\Delta x}=\left(\mathrm{b}-\mathrm{a}\right)/{\mathrm{N}}_{\mathrm{m}}$$E15

We define:

$${\mathrm{f}}_{\mathrm{i}}=\mathrm{f}\left({\mathrm{x}}_{\mathrm{i}}\right)\text{for}\mathrm{i}=0,\dots .,\mathrm{m}\text{and}:{\mathrm{x}}_{\mathrm{i}}=\mathrm{a}+\mathrm{i}.\mathrm{\Delta x}$$E16

We define the sequence:

$${\mathsf{S}}_{\mathsf{0}}=\mathsf{0}\mathsf{\text{and}}:{\mathsf{S}}_{\mathsf{i}+\mathsf{1}}={\mathrm{S}}_{\mathrm{i}}+\mathrm{\Delta x}\left(\mathrm{f}\left({\mathrm{x}}_{\mathrm{i}}\right)+\mathrm{f}\left({\mathrm{x}}_{\mathrm{i}-1}\right)\right)/2$$E17

and the sequence:

$${\mathrm{rx}}_0=0\mathrm{et}{\mathrm{rx}}_{\mathrm{i}}={\mathrm{S}}_{\mathrm{i}}/{\mathrm{S}}_{\mathrm{m}}$$E18

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

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

4.2.3.6 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 SiH₄/H₂ 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 SiH₄/H₂ 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₂), and silyl (SiH₃). The hydrogen consumption reactions SiH₄ + H → SiH₃ + H₂ and SiH₃ + H → SiH₂ + H₂ are found to play a central role in deciding the distribution of hydrogen [39].The plasma chemistry indicates that H atoms and SiH₃ radicals play an important role in the a-Si:H deposition process [40]. Experimentally, it is generally accepted that SiH₃ radicals dominate a-Si:H and μc-Si film growth from SiH₄ plasmas in the PECVD; it is the key precursor of a-Si:H deposition [41]. The proposed MCS allowed to get the ratio SiH₂/SiH₃ and mean value of densities of species. It provides information on SiH₄ dissociation and on the production of SiH₃, H, SiH₂, and Si₂H₆ and other important parameters.

The plasma in the PECVD reactor is weakly ionized. For our study, the mixture gas contains 22% of SiH₄ and 78% of H₂; 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⁸ cm⁻³. The process is considered to be stationary. We take into account electrons and eight neutral species (SiH₄, SiH₃, SiH₂, H, H₂, Si₂H₆, Si₂H₅, 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.

Symbol	Reactions	Kreac (cm^3/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 × 10^10× (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 × 10^10 × (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]

Table 1.

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

Let $K_{\mathit{\text{reac}}}^{\mathit{\text{cons}}}$ and $K_{\mathit{\text{reac}}}^{\mathit{\text{prod}}}$ be the rate constants of the consumption and the production of species A. The chemical reaction for the consumption of A is as:

$$a·A+b·B{\displaystyle \stackrel{K_{\mathit{\text{reac}}}^{\mathit{\text{cons}}}}{\to }}c·C+d·D$$

And chemical reaction for the production of A is as:

$$a\text{'}·A\text{'}+b\text{'}·B\text{'}{\displaystyle \stackrel{K_{\mathit{\text{reac}}}^{\mathit{\text{prod}}}}{\to }}c\text{'}·A+d\text{'}·D\text{'}$$

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

$$R_A=+\sum _{K_{\mathit{\text{reac}}}^{\mathit{\text{prod}}}}{K}_{\mathit{\text{reac}}\left(A^{\prime }{B}^{\prime }\right)}^{\mathit{\text{cons}}}{n}_{A\text{'}}^{a\text{'}}{n}_{B\text{'}}^{b\text{'}}-\sum _{K_{\mathit{\text{reac}}}^{\mathit{\text{cons}}}}{K}_{\mathit{\text{reac}}\left(A,B\right)}^{\mathit{\text{cons}}}{n}_A^a{n}_B^b$$E19

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

Figure 5.

Form of the simulation cell.

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:

$${Pr}_{\mathrm{sp},\mathrm{i}}=\frac{n_i}{{\sum }_{j=1}^8{n}_j}$$E20

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.

5.2.2 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:

$${\upsilon }_{\mathrm{ij}}=V_{\mathit{ij}}{n}_j<{\sigma }_{\mathit{ij}}>$$E21

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

The mean free path <λ_ι> of species i is:

$$<{\lambda }_{\mathrm{i}}>=\frac{1}{n_j<{\sigma }_{\mathit{ij}}>}$$E22

The time between two collisions τ_ij is then:

$${\tau }_{\mathit{ij}}=\frac{<{\lambda }_{\mathrm{i}}>}{V_{\mathit{ij}}}=\frac{1}{{\upsilon }_{\mathit{ij}}}$$E23

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

$$k_{\mathit{ij}}=<{\sigma }_{\mathit{ij}}\left(V_{\mathit{ij}}\right)·V_{\mathit{ij}}>$$E24

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:

$${Pr}_{\mathrm{col},\mathrm{j}}=\frac{{\nu }_{\mathit{ij}}}{{\sum }_{k=1}^9{\nu }_{\mathit{ik}}}$$E25

where $\nu$_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:

$$E_a=-k_B\mathit{Tln}\left(K_{\mathit{\text{reac}}}/{\nu }_{\mathit{ij}}\right)$$E26

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₄ 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):

$${Pr}_{\mathit{\text{reac}}}\left(i,j\right)=\frac{K_{\mathit{\text{reac}}}\left(i,j\right)}{\sum K_{\mathit{\text{reac}}}\left(i,j\right)}$$E27

where $\sum K_{\mathit{\text{reac}}}\left(i,j\right)$ 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:

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

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

Table 2.

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

For radicals (e.g., SiH₃), 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 * 10⁹; Qnp5 = 0.20296956 * 10⁹; 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 Si₂H₆, SiH, and Si₂H₅ particles reaching the surface is negligible.

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

• 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₂ and SiH₄ by R5, R1, and R2 reactions.

• The production of SiH₃ is done by R8, and then there is production of SiH₂ by R12.

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

• The second important chemical reaction in the SiH₄ dissociation is R1: SiH₄ + e → SiH₃ + 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

Table 3.

Ratios N_s,i/N_s, H₂ of particles reaching the surface compared to H₂.

Type	SiH4	SiH3	SiH2
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.

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 SiH₄/H₂ 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.