Monte Carlo Simulation of Epitaxial Growth

Abstract

A numerical Monte Carlo (MC) model is described in detail to simulate epitaxial growth. This model allows the formation of structural defects, like substitutional defects and vacancies, and desorption of adsorbed atoms on the surface. The latter feature supports the study of epitaxial growth at very high kinetic regime. The model proposed here is applied to simulate the homoepitaxial growth of Si. The results obtained fit well to the experimental reports on (0 0 1) silicon homoepitaxy. The easy implementation of a large number of microscopic processes and the three-dimensional spatial information during the film growth suggests that the model can be applied to simulate the growth of binary, ternary, or more compounds and even the growth of superlattices and heterostructures.

1. Introduction

Computer simulation has been successfully applied to the study of surface growth by molecular beam epitaxy (MBE) [1]. These simulation results play a key role in understanding and elucidating various aspects of MBE growth. These simulations provide an atomistic interpretation of the changes in characteristic reflection high energy electron diffraction (RHEED) patterns, observed in real-time, that are related to the different growth modes, i.e., from island nucleation to layer-by-layer growth mode. Additionally, by simulation, the dependency of substrate temperature and the growth rate in the vacancy formation and structural defects can be evaluated [2, 3, 4, 5, 6]. However, the majority of these computational models for simulating neglects desorption effect and, consequently, do not explain growth at high substrate temperatures, where the growth rate is decreased due to desorption of atoms adsorbed on the surface. Besides, the models often use structural approximations, in which the formation of defects is not allowed. In this work, a numerical Monte Carlo (MC) model is presented, in which desorption processes and structural defects are allowed, being possible to study limit cases, i.e., low kinetic energy regime, where structural defects are more likely, and high kinetic energy regime, where the desorption rate competes with the deposition rate.

This chapter is organized as follows: Section 2 presents a review on molecular beam epitaxy, focusing on the growth modes and the most common experimental techniques used to characterize thin epitaxial films. Section 3 presents a brief review of computational models of epitaxial growth simulation, and the model developed in this work. Section 4 presents the simulation results of silicon (Si) homoepitaxial growth on Si (0 0 1). Finally, Section 5 presents the main conclusions and new applications of the proposed model.

2. Molecular beam epitaxy

Molecular beam epitaxy (MBE) is a state-of-the-art ultra-high vacuum (UHV) thin film growth technique. The materials that will compose the films are sublimated from highly stable effusion cells, forming the molecular beam, which is deposited on an independently heated substrate. The substrate is a monocrystalline material and the impinging atoms, from the solid source, follow the substrate crystalline orientation. The physicochemical interaction mechanisms between two phases, in this case, vapor solid, with the growing solid, is called epitaxy, which has a Greek root: epi means “above” and taxis means “an ordered manner” [7].

Since the system operates in a UHV environment, the mean free path of the vapor species is much larger than the distance between the solid sources and the substrate. The interactions of the molecular beam, before colliding with the substrate surface, can be neglected. The substrate temperature is controlled independently of the effusion cells and is kept low relative to the temperatures of the cells. Therefore, the growth far from the thermodynamic equilibrium makes it possible to compensate distinct thermodynamic properties of the different materials, besides allowing the growth of high-quality monocrystalline films.

The UHV conditions in the MBE growth chamber result in an extremely low background impurity level, which allows the growth of samples with high doping control. Besides, several materials can be sublimated simultaneously and, by means of independent shutter control, multi-layer systems, with very sharp interfaces, can be obtained. Among that, a wide range of alloy composition and different doping levels can be achieved. The growth conditions are reproducible, highly stable and can vary over a wide range, which is crucial for optimizing the growth for specific materials. Additionally, compatibility with in situ analysis tools provides insights into the microscopic processes involved in the growth. The reflection high-energy electron diffraction (RHEED), for example, probes the film surface and, due to the small angle of the incident electron beam, provides information on surface morphology, in real time, during the growth [8].

In the past decades, the interest in MBE was promoted mainly by the exciting properties of semiconductor structures due to two characteristics: high control in the atomic level and reproducibility. Nowadays, after the theoretical prediction of the topological phase of the matter, the MBE technique has shown to be a promising way to obtain high-quality samples, without using counter doping to suppress free carriers due to the structural defects [9].

2.1. Growth mechanisms

Growth of MBE can take place on a substrate composed of the same material, e.g., silicon epitaxy on Si substrates [10], or on different materials, e.g., Ge epitaxy on Si [11]. The first is called homoepitaxy and the second is called heteroepitaxy. These growth conditions may lead to different growth mechanisms, due to the differences in the lattice parameters. Basically, there are three different modes of growth: Volmer-Weber, Stranski-Krastanov, and Frank-van der Merwe, which will occur depending on the experimental parameters and the lattice mismatches.

In Volmer-Weber’s mode, the interaction of the adsorbed atoms is much stronger among them than with the substrate surface, which leads to the formation of clusters or three dimensional islands. As growth proceeds, these islands expand, widening their volumes, whose height greatly exceeds the thickness of a monolayer (ML), leading to the simultaneous growth of atomic layers with rough surface, as shown in Figure 1(a) [8].

Figure 1.

Modes of growth: (a) Volmer-Weber; (b) Frank-van der Merwe; and (c) Stranski-Krastanov.

In Frank-van der Merwe’s mode, the atoms adsorbed on the surface have a stronger interaction with the surface, leading to the formation of a complete ML before another starts to grow. This layer-by-layer growth mode is often referred as a bi-dimensional growth and is shown in Figure 1(b) [8].

The last mode of growth, called Stranski-Krastanov, is characterized by mixing both previously mentioned. In this mode, the adsorbed atoms begin to grow layer-by-layer. When the critical layer thickness is reached, which value depends on the specific physical properties of each compound, the elastic energy accumulated in the growth relaxes, resulting in the formation of clusters or islands. This mode of growth is shown in Figure 1(c) [8].

Of course, all these growth modes are important, since each has a particular application. For the growth of topological insulators, such as bismuth chalcogenide compounds, the Frank-van der Merwe mode is necessary, because the aim here is to minimize structural defects during the formation of the epitaxial layers [12]. However, to manufacture low-dimensional structures, such as quantum dots, when island morphology is required, the growth parameters of the Volmer-Weber mode are determined to control the density, size, and distribution of the islands.

2.2. Characterization techniques

RHEED and atomic force microscope (AFM) are the most common experimental techniques applied to characterize the surface of epitaxial films. The first, as mentioned, provides information in situ and in real time on surface reconstruction, growth rate, and growth mode. The second allows probing the surface morphology of the film.

The RHEED equipment, basically, consists of an electron gun, in the energy range of 10–50 keV, and a phosphorescent screen. The electron beam is directed to the sample surface, at low angle (<5°). The de Broglie wavelength of electrons for this energy range is 0.17–0.06 Å, which corresponds to the interatomic distances in the crystalline lattice. Therefore, whenever the difference between the incident beam and the diffracted one is equal to a vector of the reciprocal lattice, there will be a maximum of diffraction.

The AFM technique provides high-resolution images of the surface morphology. A piezoelectric is used to scan the sample surface, and a very sharp tip is used to probe the surface. The tip deflection is measured at each point, providing spatial information in the real space of the surface. This information can be directly compared to the simulation results, in order to validate the growth models.

3. Molecular beam epitaxial models

3.1. Epitaxial growth models

The most common nucleation models used in the simulation of growth by MBE are either completely deterministic or totally stochastic. The deterministic models are based on the temporal evolution of differential equations and study microscopic parameters or stability properties of the growing surface. The deterministic models do not contain explicit spatial information of the growing surface [13, 14, 15].

Alternatively, to the analytic simulations, there are models that consider the stochastic nature of the microscopic processes. These stochastic calculations are typically implemented in the form of molecular dynamics (MD) or kinetic Monte Carlo (MC) simulations [16, 17, 18]. The MC simulation allows the easy implementation of a large number of microscopic processes. The rate of each microscopic process is obtained from first principle calculations or from experiments. These models are often in agreement with experimental works and provide spatial information of the growing surface [19].

The solid-on-solid (SOS) approach is very common in Monte Carlo simulation models [20]. In this approximation, one atom can only be accommodated on another atom, and therefore, structural defects, such as vacancies, are not allowed. Since bulk defects are not allowed in the SOS model, the growth quality can only be evaluated by the growing surface. This approach has been widely used in growth simulations of MBE due to the successful interpretations of RHEED intensity oscillations and the surface atomistic processes.

3.2. The Monte Carlo epitaxial model

This MC computational growth model was implemented using the nearest-neighbor and lattice-gas approximation. The first approximation, as mentioned, considers the interaction only with the in-plane nearest neighbor. Figure 2(a) shows an example for a cubic lattice structure, which has a maximum of four in-plane neighbors. Therefore, the potential energy of each lattice point is determined by analyzing the number of bonds between the closest atoms. This approximation has been widely used in MBE simulations [21]. The second approximation considers a fixed crystalline structure, with sites at which the atoms from the solid source can be accommodated. In this case, the possible crystalline positions can be either occupied or empty. Once a fixed crystalline structure is considered, no strain is accounted in this model. Figure 2(b) illustrates the lattice gas approximation. It is good to emphasize that these approximations do not allow studying surface reconstruction, since the atomic positions are pre-determined and the atoms can fill or not each one. Besides, these approximations allow only the study of homoepitaxy or, in some cases, heteroepitaxy of lattice-matched substrates, e.g., the growth of bismuth telluride on BaF₂ (1 1 1), where the lattice mismatch is of only 0.04% [12]. In similar cases, the lattice gas approximation can be employed in the numerical models.

Figure 2.

(a) Determination of the number of nearest neighbor to a chosen position of a cubic lattice. In this example, the atom has three in-plane bonds and one out-of-plane bond. (b) The lattice gas approximation considers a fixed empty crystalline structure, in which the atoms from the solid source are accommodated, initiating epitaxial growth.

The processes considered, in the model presented in this chapter, are shown in Figure 3. After deposition (a), surface atoms are allowed either to migrate to another position (b–e) or to desorb (g), and structural defects can be formed (f). The surface is defined as an occupied lattice position with dangling bonds. This statement is also valid for atoms near a hole, since these atoms have one or more dangling bonds due to the presence of the hole itself. This means that the neighboring atoms can migrate and fill the hole, turning on the bulk diffusion process.

Figure 3.

Allowed processes during MC simulation of thin film growth by MBE: (a) deposition of atoms, arriving from the solid sources, on the surface of a growing film; (b) surface diffusion to a more stable position; (c) atom in an unstable position, without lateral bonds, which can either desorb to the vacuum volume or diffuse to a more stable position; (d) surface diffusion or migration to a higher kink lattice; (e) possible surface diffusion forming overhang; (f) hole formed in the structure volume; and (g) desorption.

Models in which bulk defects are considered allow for the study of structural defect formation at low kinetic energy regimes, e.g., growth by MBE at low substrate temperatures. In addition, the desorption process of atoms at the surface is critical for investigating high-energy growth regimes, since desorption of adsorbed atoms (adatoms) significantly affects the growth rate and morphology of films at very high substrate temperatures.

The implementation of this model can be divided into two steps: deposition and surface analysis. The first step is responsible for the seating of the atoms on the growing surface, which, in a computational view, can be seen as a square matrix being filled with numbers. In sequence, the second step starts, and each position of the matrix is analyzed. At this point, for each matrix position, a calculation is made and, depending on the result, one of the possibilities described in Figure 3 can occur, except for the process shown in Figure 3(a).

These steps are repeated until a certain growth time, as shown in Figure 4. For each loop interaction, a time unit is added to the deposition time. By using deposition rates equivalent to the experimental ones and increasing the number of site-analysis during the second step, i.e., the number of times each position is analyzed, it is possible to use an approximately equivalent experimental time unit.

Figure 4.

Flow chart showing the two-step model. During one time unit, a number of atoms are randomly deposited on the growing surface. After deposition, all surface sites are randomly analyzed. Surface diffusion or desorption can occur during this process. The deposition time is increased by one unit, and the process is repeated until the total growth time is achieved.

3.2.1. Deposition rate

Depending on the temperature of the effusion cells, different growth rates are achieved. These rates can be experimentally determined, in situ, by measuring RHEED intensity oscillations or using a quartz crystal microbalance. In addition, these rates can be determined by directly measuring the film thickness after growth.

In the model, a certain amount of atoms is deposited on the first step of each round of the simulation. This amount depends directly on the substrate area, as given by Eq. (1),

$$N_X=L^2.D$$E1

In this equation, N_X is the number of atoms from specie X that will be deposited per second, and L is the lateral dimension of the substrate (the size of one side of the square matrix). The constant D is the reciprocal of the experimental deposition rate, given in seconds per mono layers (s/ML), i.e., D = 1/GR, where GR is the experimental growth rate in ML/s. In this model, it is possible to simulate the growth of materials composed by one or more atomic species. In that case, the first step of the model will be repeated N_tot times, with N_tot = N_X + N_Y + …, where X, Y indicate each atomic species. To fit the simulation to the experimental data, D must be carefully determined for each atomic species in the model.

3.2.2. Surface analysis

In the second step, each position of the surface is analyzed randomly. The surface control is managed by a linear dynamic data structure (implemented through a linked list), containing all spatial coordinates and the atomic species.

In a certain drawn position, the probability of surface diffusion or desorption is calculated. These probabilities are given as the product of an observation time, τ, and a hopping rate, R_{D, E}, for diffusion (D) or desorption (E), as shown in Eq. (2),

$$P_{D,E}=\tau .R_{D,E}$$E2

Theoretically, the observation time should be as low as possible, τ → 0. However, a very low observation time leads to a high computational time, since the second step is repeated until the integrated observation time is equal to one unit, as shown in Eq. (3). In practice, a good value is τ ≈ 0.01, which means that the second step of the model is repeated 100 times,

$$\sum _{i=1}{\tau }_i=1.$$E3

The hopping rate is determined by an Arrhenius equation, as shown in Eq. (4),

$$R_{D,E}=R_0.e^{-\frac{E_{D,E}}{k_{B}T}}$$E4

In this equation, E_D,E is the energy for diffusion (D) or desorption (E), k_B is the Boltzmann constant, and T is the substrate temperature. The typical vibration frequency of the atom,R₀, is function of Planck constant, h, and is given by Eq. (5), typically around 10¹³ Hz,

$$R_0=\frac{2k_{B}T}{h}\text{.}$$E5

The energies for diffusion or desorption can be obtained from ab initio calculations or by experimental results. These values can be fitted to the experimental data to reproduce the experimental results.

The energy depends on the number of bonds that each atom possesses. Equation (6) presents the energy required for an atom to diffuse (E_D) or desorb (E_E),

$$E_{D,E}=m{E}_{D0,E0}+n{E}_{DL,\mathit{EL}}\text{,}$$E6

where m is the number of out-of-plane bonds, n is the number of in-plane bonds (see Figure 2(a)), and the sub-indices D0 and E0 indicate the energy required to diffuse or desorb, respectively, for one out-of-plane bond.

The sub-indices DL and EL indicate the energy required to diffuse or desorb, respectively, for one in-plane bond. In practice, in a cubic lattice, m = 0 or 1 and n = 0, 1, 2, 3 or 4. When an atom is located on a vacancy, which is an unoccupied site in the lattice, m = 0.

3.2.3. Probabilities

At each round, in the second step, all surface positions are analyzed once. For each position, the number of in-plane and out-of-plane bonds is computed, and then, the probabilities for diffusion and desorption are calculated. A random number is generated, and if it reaches one of these events, the atom will migrate to another available position or desorb. Figure 5(a)shows a draw number, where none of the events was reached, and Figure 5(b) shows the activation of a surface diffusion event.

Figure 5.

(a) Example of a draw number where none of the events is reached and (b) when a surface diffusion event is achieved.

When a diffusion event is reached, the atom will migrate to any random neighboring position in the lattice. For a desorption event, the atom is simply removed from the surface. Since the MBE equipment operates under UHV conditions, the mean free path inside the growth chamber is very high, and these atoms can stick to the MBE wall chamber without interacting with the molecular beam.

For very high substrate temperatures, some lattice positions can be considered unstable, since the probability of diffusion and/or desorption can surpass one unit. Whenever the sum of the probabilities surpasses one unit, as shown in Eq. (7), a “right event” is achieved,

$$P_D+P_E>1\text{.}$$E7

This situation tries to describe limiting conditions, where the growing rate is overlapped by the desorption rate. Since this event occurs at high substrate temperatures, the adatoms have a very high surface mobility. To accomplish this condition, each time a “right event” is reached all available surface positions are analyzed. If one or more of these positions offers enough bonds to avoid another “right event,” the atom migrates randomly to one of these positions. Otherwise, if none of the surface positions is stable enough, i.e., if Eq. (6) is true, the atom desorbs to the vacuum.

To avoid border effects, a periodic boundary condition was implemented. Atoms located at the borders of the substrate have only two or three in-plane bonds. By increasing the substrate temperature, these positions can become more unstable. The toroid considers a closed substrate, such as: x(1) ↔ x(L) and y(1) ↔ y(L). In this sense, if an equivalent position is occupied, x(1) ↔ x(L), the atom located at x(1) gets a bond.

3.2.4. Random numbers

To obtain the stochastic nature of the microscopic processes, the Mersenne Twister random number generator was employed, which is a 623 dimensionally equidistributed uniform pseudorandom number generator with an accuracy of up to 32 bits [22].

3.2.5. Quantifying results

Surface roughness is estimated based on the exposed surface, without computing any structural defects. The simulated film roughness is calculated by:

$$\sigma =\frac{1}{L^2}\sum _{x=0}^L\sum _{y=0}^L\sqrt[2]{(h_{(x,y)}-\overline{h})2}.$$E8

The term h_(x,y) represents the highest position occupied by an atom, and $\overline{h}$ is the average thickness of the entire film. The equation for roughness estimation resembles the R_sq coefficient, obtained experimentally from atomic force microscopy (AFM), scanning tunneling microscopy (STM), or X-ray reflectivity measurements, which allows a direct comparison between simulation and experimental results reported in the literature (Section 2.2).

RHEED patterns are very sensitive to the surface roughness and morphology of the outer layers. The growth rate can be estimated in situ during MBE growth through the RHEED intensity oscillations, which is dynamically calculated by:

$$I\left(t\right)=(\sum _n{(-1)}^n.S_n\left(t\right))2.$$E9

where the term S_n(t) is the exposed cover of the nth layer at time t. This equation is slightly different from the equations used in the SOS models, since it calculates the reflected intensity using only surface atoms exposed to the beam [23].

4. Simulation results

4.1. Silicon (0 0 1)

For the simulation of a (0 0 1) Si homoepitaxy, the energy value for diffusion to out-of-plane bond is E_D0 = 1 eV and the in-plane energy is E_DL = 0.5 eV. The energy barrier for the evaporation process was determined by fitting the growth rate curves as a function of substrate temperature with the calculated growth rates. The curves were calculated using the vapor pressure of silicon as a function of temperature and the Knudsen equation [8]. The out-of-plane energy barrier found for evaporation is E_E0 = 3.8 eV and the in-plane energy is E_EL = 0.3 eV. The out-of-plane barrier energy determines the inflection point of the growth rate curve, and the in-plane energy determines directly the derivative of the curve, which is the rate of change as a function of substrate temperature.

During growth simulation, the exposed surface, surface roughness, RHEED intensity, structural defects density, growth rate, and the partial structure, containing all atoms coordinates, were continuously recorded as a function of time. The substrate temperature was investigated in a range varying from room temperature up to the silicon melting point. The external parameters of the simulation are substrate temperature, deposition rate, lattice size, and deposition time.

Figure 6 presents surface roughness in monolayers (ML) for four different substrate temperatures (T_SUB). The deposition rate was kept constant at 0.01 ML/s during 1000 s of deposition, which results in films thicknesses with approximately 10 ML. At room temperature (300 K), the atoms on surface do not have sufficient energy to migrate leading to a limited surface diffusion, which favors the formation of structural defects, like vacancies and overhangs. As the substrate temperature increases to 600 K, the diffusion hopping rate raises and atoms situated on unstable positions, with a few bonds, may migrate to more stable positions on the lattice. This process decreases the surface roughness, as observed in Figure 6, since the innermost layers tend to be fulfilled. The growing surface at this condition is dominated by coalescence of small clusters. For substrate temperature of 800 K, a layer-by-layer growth mode is reached. This condition is evidenced by very clear oscillations on the surface roughness profile, with a period of approximately 100 s. The oscillation ends in a minimum, indicating that the 10th layer is practically completed. Increasing the substrate temperature to 1200 K, the surface roughness is raised to a value around 1 ML, which is a consequence of film thickness inhomogeneity, due to high surface diffusion rates.

Figure 6.

Surface roughness profiles for films grown on different substrate temperatures (T_SUB). Due to limited surface diffusion for films grown at room temperature, high roughness is obtained. Raising the substrate temperature, a high surface mobility is achieved, leading to a step-by-step growth mode at 800 K. For even higher substrate temperatures, atoms are allowed to diffuse with elevated rate and the reevaporation becomes significant, enhancing surface roughness.

The results of RHEED intensity oscillations dynamically calculated for a set of growth conditions are displayed in Figure 7. At room temperature, the calculated intensity does not exhibit any oscillations, which indicates that the intensity of the interference pattern is suppressed by the rough surface. The growth at T_SUB between 600 and 800 K exhibits well-defined intensity oscillations with a period around 100 s, indicating a growth rate of 1 ML/s. At 600 K, the oscillations are less intense and less symmetric than at 800 K, indicating the growth of more than one layer at the same time. For substrate temperatures around 800 K a layer-by-layer growth mode is achieved, evidenced by the symmetric RHEED intensity oscillation and confirmed by the calculated surface roughness (Figure 6). At 1200 K, the intensity oscillations disappear again. The layer-by-layer mode is lost due to higher surface atoms mobility, which leads to rougher surfaces. These results agree very well to experimental work of (0 0 1) Si homoepitaxy [24]. In this report, the RHEED intensity oscillations are weak and vanish after 2 min of deposition at room temperature and at 1270 K, whereas they are intense, well defined, and longstanding for substrate temperatures between 600 K and 900 K. Surface coverage of the first 10 layers was calculated and is presented in Figure 8(a) as a function of substrate temperature.

Figure 7.

Normalized RHEED intensity oscillations dynamically calculated for several growth substrate temperatures (T_SUB). The intensity oscillations are suppressed at room temperature and at 1200 K, indicating rough surfaces at these conditions. For intermediated substrate temperatures, from 600 to 800 K, clearly defined oscillations are observed.

Figure 8.

(a) Surface coverage as a function of deposition time of the first ten atomic layers for substrate temperatures from 300 K to 1200 K. (b) Exposed surface for a L = 100 lattice positions and four substrate temperatures. Increasing substrate temperature enhances atoms mobility, lowering surface roughness by widening terraces on the surface (T_SUB from 300 to 700 K). At T_SUB = 800 K, a transition from cluster to layer-by-layer growth is achieved.

At lower temperatures, more than one atomic layer is filled at the same time, due to low surface diffusion. This is evidenced in the lower panel of Figure 8, where the 10th layer starts to be filled around 400 s, whereas the other 9 layers are still been filled. At T_SUB = 600 K, the coverage chart shows that simultaneous growth occurs for a maximum of two layers, whereas, at 800 K, the chart evidences the growth of only one layer per time period. Raising T_SUB to 1200 K, the high diffusion hopping rate recovers simultaneous layer growth and generates noisy lines on the surface coverage profile. Besides that reevaporation process becomes significant, contributing to an increase of spiked profile due to an abrupt change on the atomic layers coverage.

The exposed surface is presented in Figure 8(b) for growth temperatures between 300 K and 800 K.

From room temperature to T_SUB = 500 K, the images present a rough surface with cluster of Si atoms close to each other. At 700 K, the surface becomes much more flat, exhibiting large terraces with monolayer islands and a few void lakes. This growth condition is mostly controlled by islands coalescence and is on transition to the layer-by-layer growth regime. Raising the substrate temperature to 800 K, a plain layer-by-layer growth mode is observed, where the lowermost atomic layers are completely fulfilled by a step-flow mechanism before a new layer start to be formed. These results are in agreement with experimental STM images captured in an ultra-high vacuum MBE system right after the (0 0 1) Si homoepitaxy growth [25].

The dependence of the defects density formed during film growth as a function of the substrate temperature is presented in Figure 9 for different deposition rates. This graph shows that, for each deposition rate, there is a threshold temperature where the defect density starts to decay exponentially. This temperature is higher for higher deposition rates. As the deposition rate increases, higher diffusion rates, i.e., higher temperatures are needed to accommodate the impinging atoms without the formation of holes. The local minima observed in the curves of Figure 9 correspond to the variation of the diffusion rate due to the average number of bonds of the surface atoms.

Figure 9.

Density of holes formed in the simulated films as a function of substrate temperature for deposition rate varying from 0.01 ML/s up to 2.56 ML/s.

At lower substrate temperatures regimes, the hopping rate is decreased, inhibiting the atoms mobility. These growth conditions are favorable to structural defects, which leads to rough surfaces. For higher temperatures, in an energetically unstable condition with many dangling bonds, the hopping rate for diffusion raises. This condition favors atoms located in unstable positions to search other points in the lattice with lower energy, decreasing the number of structural defects and the surface roughness. With higher substrate temperatures, the atoms located in the surface are allowed to migrate in a short period of time. This condition allows atoms to migrate to a more stable position, like on kinks, giving rise to a step-flow growth mode with flatter surfaces and less structural defects. For even higher substrate temperatures, the hopping rate for diffusion is very high, even for positions with many bonds. In this situation, several diffusion movements are allowed, producing inhomogeneous thicknesses with rougher surfaces. The hopping rate for reevaporation begins to be significant, lowering the growth rate. If substrate temperature continues to increase, the reevaporation rate becomes more and more significant, until growth rate vanishes.

4.2. Applications

This model can be applied to study low and high kinetic growth regimes, since this numerical model of MBE growth considers bulk defects and desorption process and the growth of multiple elements, such as binary and ternary alloys. Besides, it is possible to study thermal treatment effects on crystalline structures, sputtering processes of surfaces, and even the growth on substrates with cleavage steps.

One topic example, which has attracted a lot of attention from experimental scientists, is topological insulator materials. The most common archetypes of this new electronic state of matter are the bismuth chalcogenides compounds, such as Bi₂Te₃ and Bi₂Se₃. This class of matter has a nontrivial topological order, resulting in a metallic behavior attributed to the surface states but insulating behavior in the bulk. However, the presence of spontaneous structural defects in these compounds leads to electric conduction in the bulk, which overlaps with surface states. Of course, in order to achieve practical applications, it is fundamental to have insulating samples in the bulk. The model presented here can be applied to understand the growth dynamics of these compounds, elucidating the mechanisms of structural defect formation.

5. Conclusion

This chapter presented firstly a brief description of the molecular beam epitaxy and the main growth modes achieved using this technique as well as of the most common experimental tools to probe thin film surfaces. The second part described a numerical model based on Monte Carlo method to simulate epitaxial growth.

The model proposed here is applied to simulate the molecular beam epitaxial growth of Si on (0 0 1) Si substrate. The results obtained fit well to the experimental reports on (0 0 1) silicon homoepitaxy. At low substrate temperatures, films exhibit a high density of structural defects and a very rough surface. On the other hand, for high temperatures, large surface terraces are achieved with a low density of bulk defects. At even higher substrate temperatures, the surface roughness increases due to the high mobility of the atoms at the surface, which could lead to non-homogeneity of the film thickness. By further increasing the substrate temperature, the atoms have sufficient energy to desorb, resulting in a reduction of the growth rate.

The large number of microscopic processes and the three-dimensional spatial information during the film growth allow the model to simulate the growth of binary, ternary, or more compounds and even the growth of superlattices and heterostructures. In addition, the model allows an easy implementation to study thermal treatments in crystalline structures and surface sputtering processes.