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.
Keywords
- Monte Carlo simulation
- molecular beam epitaxy
- epitaxial growth
- lattice-matched substrates
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:
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 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].
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
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 BaF2 (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.
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.
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.
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 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),
In this equation,
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,
Theoretically, the observation time should be as low as possible,
The hopping rate is determined by an Arrhenius equation, as shown in Eq. (4),
In this equation,
The energies for diffusion or desorption can be obtained from
The energy depends on the number of bonds that each atom possesses. Equation (6) presents the energy required for an atom to diffuse (
where
The sub-indices
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.
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,
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:
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:
The term
RHEED patterns are very sensitive to the surface roughness and morphology of the outer layers. The growth rate can be estimated
where the term
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
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 (
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
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
The exposed surface is presented in Figure 8(b) for growth temperatures between 300 K and 800 K.
From room temperature to
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.
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 Bi2Te3 and Bi2Se3. 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.
References
- 1.
Madhukar A, Ghaisas SV. The nature of molecular beam epitaxial growth examined via computer simulations. Critical Reviews in Solid State and Materials Sciences. 1988; 14 :1-130. DOI: 10.1080/01611598808241266 - 2.
Family F, Lam PM. Renormalization-group analysis and simulational studies of groove instabilitiy in surface growth. Physica A: Statistical, Mechanics and Its Applications. 1994; 205 :272-283. DOI: 10.1016/0378-4371(94)90504-5 - 3.
Levi AC, Kotrla M. Theory and simulation of crystal growth. Journal of Physics: Condensed Matter. 1997; 9 :299-344. DOI: 10.1088/0953-8984/9/2/001 - 4.
Das Sarma S, Tamborenea P. A new universality class for kinetic growth: One-dimensional molecular-beam epitaxy. Physical Review Letters. 1991; 66 :325-328. DOI: 10.1103/PhysRevLett.66.325 - 5.
Das Sarma S, Lanczycki C, Kotlyar R, Ghaisas S. Scale invariance and dynamical correlations in growth models of molecular beam epitaxy. Physical Review E. 1996; 53 :359-388. DOI: 10.1103/PhysRevE.53.359 - 6.
Landau DP, Pal S. Monte Carlo simulation of simple models for thin film growth by MBE. Thin Solid Films. 1996; 272 :184-194. DOI: 10.1016/0040-6090(95)06945-3 - 7.
Bauer G, Springholz G. Molecular beam epitaxy—Aspects and applications. Vacuum. 1992; 43 :357-365. DOI: 10.1016/0042-207X(92)90038-X - 8.
Herman MA, Sitter H. Molecular Beam Epitaxy: Fundamentas and Current Status. 2nd ed. Berlin: Springer; 1988. p. 382. DOI: 10.1007/978-3-642-97098-6 - 9.
Ando Y. Topological insulator materials. Journal of the Physical Society of Japan. 2013; 82 :102001-1-102001-32. DOI: 10.7566/JPSJ.82.102001 - 10.
Schelling C, Springholz G, Schäffler F. New kinetic growth instabilities in Si(001) homoepitaxy. Thin Solid Films. 2000; 369 :1-4. DOI: 10.1016/S0040-6090(00)00823-3 - 11.
Eaglesham DJ, Cerullo M. Low‐temperature growth of Ge on Si(100). Applied Physics Letters. 1991; 58 :2276-2278. DOI: 10.1063/1.104898 - 12.
Fornari CI, Rappl PHO, Morelhão SL, Abramof E. Structural properties of Bi2Te3 topological insulator thin films grown by molecular beam epitaxy on (111) BaF2 substrates. Journal of Applied Physics. 2016; 119 :156303-1-156303-9. DOI: 10.1063/1.4947266 - 13.
Ratsch C, Venables JA. Nucleation theory and the early stages of thin film growth. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films. 2003; 21 :s96-s109. DOI: 10.1116/1.1600454 - 14.
Kariotis R, Lagally MG. Rate equation modelling of epitaxial growth. Surface Science. 1989; 216 :557-578. DOI: 10.1016/0039-6028(89)90395-6 - 15.
Aumann CE, Kariotis R, Lagally MG. Rate equation modeling of interface width. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films. 1989; 7 :2180-2185. DOI: 10.1116/1.575953 - 16.
Liang YY, Yoon SF, Fitzgerald EA. Kinetic Monte Carlo simulation of quantum dot growth on stepped substrates. Journal of Physics D: Applied Physics. 2013; 46 :495102. DOI: 10.1088/0022-3727/46/49/495102 - 17.
Günther V, Mauß F, Klauer C, Schlawitschek C. Kinetic Monte Carlo simulation of the epitaxial growth of Si(1 0 0). Physica Status Solidi (c). 2012; 9 :1955-1962. DOI: 10.1002/pssc.201200340 - 18.
Günther V, Mauß F. Si(100)2×1 epitaxy: A kinetic Monte Carlo simulation of the surface growth. Physics Procedia. 2013; 40 :56-64. DOI: 10.1016/j.phpro.2012.12.008 - 19.
Stumpf R, Scheffler M. Ab initio calculations of energies and self-diffusion on flat and stepped surfaces of Al and their implications on crystal growth. Physical Review B, Condensed Matter. 1996; 53 :4958-4973. DOI: 10.1103/PhysRevB.53.4958 - 20.
Weeks JD, Gilmer GH, Jackson KA. Analytical theory of crystal growth. The Journal of Chemical Physics. 1976; 65 :712-720. DOI: 10.1063/1.433086 - 21.
Marmorkos IK, Das Sarma S. Kinetic simulation of molecular beam epitaxial growth dynamics. Surface Science. 1990; 237 :L411-L416. DOI: 10.1016/0039-6028(90)90511-6 - 22.
Matsumoto M, Nishimura T. Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactionson Modeling and Computer Simulation. 1998; 8 :3-30. DOI: 10.1145/272991.272995 - 23.
Marmorkos IK, Das Sarma S. Atomistic numerical study of molecular-beam-epitaxial growth kinetics. Physical Review B. 1992; 45 :11262-11272. DOI: 10.1103/PhysRevB.45.11262 - 24.
Sakamoto T, Kawai NJ, Nakagawa T, Ohta K, Kojima T, Hashiguchi G. Rheed intensity oscillations during silicon MBE growth. Surface Science. 1986; 174 :651-657. DOI:10.1016/0039-6028(86)90487-5 - 25.
Jernigan GG, Thompson PE. Temperature dependence of atomic scale morphology in Si homoepitaxy between 350 and 800 °C on Si (100) by molecular beam epitaxy. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films. 2001; 19 :2307-2311. DOI: 10.1116/1.1384559