- Email: [email protected]
The Monte Carlo simulation of epitaxial growth of hexagonal GaN
Surface Science 432 (1999) L617–L620 www.elsevier.nl/locate/susc
Surface Science Letters
The Monte Carlo simulation of epitaxial growth of hexagonal GaN Huibing Mao, Songhoon Lee, Seong-Ju Park * Department of Materials Science and Engineering and Center for Electronic Materials Research, Kwangju Institute of Science and Technology, Kwangju 500-712, South Korea Received 20 July 1998; accepted for publication 6 May 1999
Abstract The growth mechanism of GaN was studied by Monte Carlo simulation. The simulation was carried out on the hexagonal structure of GaN. The simulation of the growth process on the vicinal surface shows that epitaxial growth is in a step-flow mode at the simulation temperatures above the critical temperature T , otherwise it is in a twoC dimensional mode. The upper limit of the energy barrier to hopping is roughly estimated to be 1.43±0.25 eV which is in a good agreement with the reported experimental value. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Density of steps; Epitaxial growth; GaN; Hexagonal structure; Island density; Monte Carlo simulation
Gallium nitride (GaN ) has been considered to be one of the most promising materials for the development of optical devices in the region of blue to UV waveband. In the past several years much progress has been made both for the materials and the devices [1–5]. Metalorganic chemical vapor deposition and molecular beam epitaxy (MBE ) have been the most common techniques for the growth of GaN and related semiconductors. However, despite many efforts made in recent years, many fundamental problems related to a growth mechanism still remain unsolved. The kinetics of GaN formation, for example, is still not well understood. Up to now, Monte Carlo simulation has been * Corresponding author. Fax: +82-62-970-2309. E-mail address: [email protected] (S.-J. Park)
widely used to study the growth mechanism of GaAs and its related semiconductors [6,7]. Monte Carlo simulation provides a simple but effective method to study the growth mechanism. For GaAs and related materials, a simple cubic lattice structure has been successfully employed in the simulation. However, the simple cubic lattice model is not suitable for a hexagonal GaN and related semiconductors which have wurtzite structures. In this study we employed a hexagonal lattice structure to model the growth processes of GaN and related semiconductors using Monte Carlo simulation. There are six adjacent sites for every adatom in the hexagonal structure whereas there are four adjacent sites for that in the simple cubic structure, which is expected to make a marked difference in the surface kinetics. For example, an adatom in a hexagonal substrate structure has a
0039-6028/99/$ – see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S0 0 39 - 6 0 28 ( 99 ) 0 06 3 9 -1
SU RF A LE CE TT S ER CIE S N CE
L618
H. Mao et al. / Surface Science 432 (1999) L617–L620
higher hopping rate than that in a cubic substrate structure under the same energy barrier to hopping. We will discuss the effect of the structural difference on the growth processes in detail. The previous studies on the growth showed that the growth mechanism of GaN is much more complex than that of other semiconductors, such as GaAs etc. It has been reported that the GaN growth is the only process which is observed at a high V/III flux ratio. However, three processes are typically observed at a low V/III flux ratio: GaN growth, Ga desorption and Ga accumulation [6 ]. In the transition region, two processes are observed: GaN growth and Ga desorption. In this study we will focus only on GaN growth process which is expected under a typical growth condition (a high V/III flux ratio), where the kinetics of Ga adatoms is the main process for epitaxial growth of GaN. In order to model the kinetics of GaN growth during the epitaxial growth, we employed a solidon-solid model [7,8]. Because the epitaxial growth is carried out under nitrogen-rich conditions, the cation kinetics, rather than the Ga–N reaction kinetics and its related one, is the rate-limiting step in the growth of GaN. The main difference between the present structure we adopted to simulate the growth of GaN and the previous one used for the growth of GaAs lies in the fact that every adatom in the hexagonal GaN has six nearest neighbors whereas that in the cubic GaAs has four nearest neighbors, as shown in Fig. 1 [9]. Therefore every
Fig. 1. A sketch of the hexagonal structure in Monte Carlo simulation, in which the six arrows indicate the six adjacent sites of an adatom.
adatom in the hexagonal GaN can hop randomly to one of its six adjacent sites with a hopping rate D, where D is the hopping rate and is determined by an Arrhenius expression: D=D exp(−E /k T ). (1) 0 r B In this expression, T is a substrate temperature, E is the energy barrier to hopping, and k is r B Boltzmann’s constant. With the assumption that each adatom has the same frequency and can be treated as a two-dimensional harmonic oscillator, a simple and reasonable approximation of D at 0 high temperatures is obtained from a straightforward application of the classical equipartition theorem: D =2k T/h, (2) 0 B where h is Plank’s constant. Because the exact value of E is not available, we adopt the hopping r rate D as a main parameter. In fact, it will be more convenient to compare the results for a hexagonal structure with those for a cubic structure to examine the parameter D. Furthermore, the energy barrier to hopping can also be estimated in the present study. To simplify our model based on the solid-onsolid model, it is also assumed that the nucleation is irreversible, which means that the smallest stable nucleus is a dimer. The statistics were obtained from ten independent Monte Carlo simulation runs on 200×200 lattices. The simulations were carried out on two kinds of substrates. The first one was a normal (0001) substrate with a periodic boundary condition. The second one was a vicinal surface which is misoriented by 2° to the [101: 0] direction from the standard [0001] direction. There are eight steps on this simulated vicinal surface and the step length is 25 lattice units. These data were obtained from the structure parameters of the wurtzite GaN. It was also assumed to have a periodic boundary in the vicinal surface in the simulation: if an adatom cross over the upper boundary in the [101: 0] direction, it will reappear in the lower boundary, and it will be the same for a reverse case. In its perpendicular direction, we used a normal periodic boundary condition. The density of steps is a major parameter which is essential to describe the growth processes. The growth mode can be determined from information
CE N IE SC RS CE E A TT RF LE SU
H. Mao et al. / Surface Science 432 (1999) L617–L620
on the density of steps. Especially, if it is compared with the reflection high-energy electron diffraction (RHEED) intensity, some important parameters, such as the energy barrier to hopping, can be determined from the variation of the density of steps. In order to compare our simulation results with some available experimental data [10], we estimated the growth processes on a vicinal surface. In the following simulation shown in Fig. 2, we take the deposition rate to be 1.66 ML per second for the comparison with the experimental results. Two distinct growth modes are noticeable in Fig. 2, depending on the hopping rates. For the hopping rate in the range of 104 to 105 s−1, the density of steps clearly shows an oscillatory behavior, implying that the GaN growth is in a two-dimensional growth mode. However, the oscillatory behaviour in the density of steps with the growth time gradually disappears with the increase in the hopping rate. For the hopping rate D=9×106 s−1; however, there is no oscillation in the density of steps, which implies that the growth mode has been changed to a step flow mode on the vicinal surface. In the RHEED intensity measurements of the GaN growth by MBE, Grandjean et al. found that the growth mode changed from the two-dimensional mode to the step flow mode at a critical temperature T =800°C and at a growth rate of C 1.66ML s−1 [10]. In their experiments the substrate misorientation was less than 2°. With a critical temperature T =800°C and a hopping rate C
Fig. 2. The density of steps on the vicinal surface for different hopping rates are shown. The vicinal surface is misoriented by 2° to the [101: 0] direction from the standard [0001] direction.
L619
D=9×106 s−1, we could roughly estimate the upper limit of the energy barrier to hopping to be 1.43 eV using the Arrhenius expression. There are several factors that affect the precision of the estimated energy barrier to hopping E : the prer factor D , the substrate misorientation and the 0 critical temperature T . The error bar of the energy C barrier to hopping E can be estimated using Eq. r (1). Among the factors mentioned, the prefactor D is the major factor leading to an error. The 0 error due to the prefactor is determined by the substrate temperature and Eq. (2), and it may lead to an uncertainty of 0.2 eV. The uncertainty in the critical temperature is less than 10°C. It only introduces an uncertainty of less than 1%. An uncertainty in the other factors is about the same. Therefore the total uncertainty in the energy barrier to hopping E is less than 0.25 eV. The r upper limit of the energy barrier to hopping is estimated to be 1.43±0.25 eV. This estimated result is very close to an experimental value of 1.45±0.25 eV observed on a Ga-terminated GaN surface by Liu et al. [11]. The activation energy of the surface migration, 1.45±0.25 eV was determined from the recovery behavior of the RHEED intensity[11]. Our model can be further refined. For example, a thermal desorption of a nucleus can be included in the model. If the thermal desorption of a nucleus is considered, the energy barrier to hopping must include two energy terms: a surface contribution E and a contribution from S the nearest neighbor along the surface E . The N values E and E can be regarded as the effective S N energy barriers to hopping. Those values incorporate, in an average way, the fast processes as well as other factors such as the surface reconstruction. According to our previous result, two parameters, E and E /E which will be included in the refined S S N model, are primarily responsible for the temperature of the substrate and the effect of lateral interactions at a given temperature, respectively [9]. It is expected that the values of E and E are S N strongly related to the growth conditions. Therefore our model based upon Monte Carlo simulation and the solid-on-solid model is a useful method to determine the energy barrier to hopping if accurate experimental results are available. In conclusion, we have employed the Monte
SU RF A LE CE TT S ER CIE S N CE
L620
H. Mao et al. / Surface Science 432 (1999) L617–L620
Carlo method and the solid-on-solid model to simulate an epitaxial growth mechanism of the hexagonal structure of GaN semiconductor. The adatoms on the hexagonal structure have more probability to hop freely than those on the cubic structure, resulting in the lower density of islands. The simulation of the growth process on the vicinal surface shows that an epitaxial growth is in a stepflow mode at simulation temperatures above the critical temperature T ; otherwise it is in a twoC dimensional mode. The upper limit of the energy barrier to hopping is estimated to be 1.43±0.25 eV, which is in a good agreement with the reported experimental result.
Acknowledgement This work was supported in part by Korea Science and Engineering Foundation ( KOSEF ).
References [1] S. Nakamura, M. Senoh, N. Iwasa, S. Nagahama, Jpn. J. Appl. Phys. 34 (1995) L797. [2] I. Akasaki, H. Amano, S. Sota, H. Sakai, T. Tanaka, M. Koike, Jpn. J. Appl. Phys. 34 (1995) L1517. [3] S. Nakamura, M. Senoh, S. Nagahama, N. Iwasa, T. Yamada, T. Matsushita, H. Kioku, Y. Sugimoto, Jpn. J. Appl. Phys. 35 (1996) L74. [4] B.W. Lim, Q.C. Chen, J.Y. Yang, M.A. Khan, Appl. Phys. Lett. 68 (1996) 3761. [5] M.A. Khan, Q. Chen, C.J. Sun, J.W. Yang, M.S. Shur, H. Park, Appl. Phys. Lett. 68 (1996) 514. [6 ] J.R. Jenny, R. Kaspi, K.R. Evans, J. Cryst. Growth 175–176 (1997) 89. [7] S. Clarke, D.D. Vvedensky, Phys. Rev. Lett. 58 (1987) 2235. [8] J.D. Weeks, G.H. Gilmer, Adv. Chem. Phys. 40 (1979) 157. [9] H.B. Mao, W. Lu, S.C. Shen, J. Cryst. Growth 151 (1995) 31. [10] N. Gradjean, J. Massies, Appl. Phys. Lett. 71 (1997) 1816. [11] H. Liu, J.G. Kim, M.H. Ludwig, P.M. Park, Appl. Phys. Lett. 71 (1997) 347.
Surface Science Letters
The Monte Carlo simulation of epitaxial growth of hexagonal GaN Huibing Mao, Songhoon Lee, Seong-Ju Park * Department of Materials Science and Engineering and Center for Electronic Materials Research, Kwangju Institute of Science and Technology, Kwangju 500-712, South Korea Received 20 July 1998; accepted for publication 6 May 1999
Abstract The growth mechanism of GaN was studied by Monte Carlo simulation. The simulation was carried out on the hexagonal structure of GaN. The simulation of the growth process on the vicinal surface shows that epitaxial growth is in a step-flow mode at the simulation temperatures above the critical temperature T , otherwise it is in a twoC dimensional mode. The upper limit of the energy barrier to hopping is roughly estimated to be 1.43±0.25 eV which is in a good agreement with the reported experimental value. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Density of steps; Epitaxial growth; GaN; Hexagonal structure; Island density; Monte Carlo simulation
Gallium nitride (GaN ) has been considered to be one of the most promising materials for the development of optical devices in the region of blue to UV waveband. In the past several years much progress has been made both for the materials and the devices [1–5]. Metalorganic chemical vapor deposition and molecular beam epitaxy (MBE ) have been the most common techniques for the growth of GaN and related semiconductors. However, despite many efforts made in recent years, many fundamental problems related to a growth mechanism still remain unsolved. The kinetics of GaN formation, for example, is still not well understood. Up to now, Monte Carlo simulation has been * Corresponding author. Fax: +82-62-970-2309. E-mail address: [email protected] (S.-J. Park)
widely used to study the growth mechanism of GaAs and its related semiconductors [6,7]. Monte Carlo simulation provides a simple but effective method to study the growth mechanism. For GaAs and related materials, a simple cubic lattice structure has been successfully employed in the simulation. However, the simple cubic lattice model is not suitable for a hexagonal GaN and related semiconductors which have wurtzite structures. In this study we employed a hexagonal lattice structure to model the growth processes of GaN and related semiconductors using Monte Carlo simulation. There are six adjacent sites for every adatom in the hexagonal structure whereas there are four adjacent sites for that in the simple cubic structure, which is expected to make a marked difference in the surface kinetics. For example, an adatom in a hexagonal substrate structure has a
0039-6028/99/$ – see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S0 0 39 - 6 0 28 ( 99 ) 0 06 3 9 -1
SU RF A LE CE TT S ER CIE S N CE
L618
H. Mao et al. / Surface Science 432 (1999) L617–L620
higher hopping rate than that in a cubic substrate structure under the same energy barrier to hopping. We will discuss the effect of the structural difference on the growth processes in detail. The previous studies on the growth showed that the growth mechanism of GaN is much more complex than that of other semiconductors, such as GaAs etc. It has been reported that the GaN growth is the only process which is observed at a high V/III flux ratio. However, three processes are typically observed at a low V/III flux ratio: GaN growth, Ga desorption and Ga accumulation [6 ]. In the transition region, two processes are observed: GaN growth and Ga desorption. In this study we will focus only on GaN growth process which is expected under a typical growth condition (a high V/III flux ratio), where the kinetics of Ga adatoms is the main process for epitaxial growth of GaN. In order to model the kinetics of GaN growth during the epitaxial growth, we employed a solidon-solid model [7,8]. Because the epitaxial growth is carried out under nitrogen-rich conditions, the cation kinetics, rather than the Ga–N reaction kinetics and its related one, is the rate-limiting step in the growth of GaN. The main difference between the present structure we adopted to simulate the growth of GaN and the previous one used for the growth of GaAs lies in the fact that every adatom in the hexagonal GaN has six nearest neighbors whereas that in the cubic GaAs has four nearest neighbors, as shown in Fig. 1 [9]. Therefore every
Fig. 1. A sketch of the hexagonal structure in Monte Carlo simulation, in which the six arrows indicate the six adjacent sites of an adatom.
adatom in the hexagonal GaN can hop randomly to one of its six adjacent sites with a hopping rate D, where D is the hopping rate and is determined by an Arrhenius expression: D=D exp(−E /k T ). (1) 0 r B In this expression, T is a substrate temperature, E is the energy barrier to hopping, and k is r B Boltzmann’s constant. With the assumption that each adatom has the same frequency and can be treated as a two-dimensional harmonic oscillator, a simple and reasonable approximation of D at 0 high temperatures is obtained from a straightforward application of the classical equipartition theorem: D =2k T/h, (2) 0 B where h is Plank’s constant. Because the exact value of E is not available, we adopt the hopping r rate D as a main parameter. In fact, it will be more convenient to compare the results for a hexagonal structure with those for a cubic structure to examine the parameter D. Furthermore, the energy barrier to hopping can also be estimated in the present study. To simplify our model based on the solid-onsolid model, it is also assumed that the nucleation is irreversible, which means that the smallest stable nucleus is a dimer. The statistics were obtained from ten independent Monte Carlo simulation runs on 200×200 lattices. The simulations were carried out on two kinds of substrates. The first one was a normal (0001) substrate with a periodic boundary condition. The second one was a vicinal surface which is misoriented by 2° to the [101: 0] direction from the standard [0001] direction. There are eight steps on this simulated vicinal surface and the step length is 25 lattice units. These data were obtained from the structure parameters of the wurtzite GaN. It was also assumed to have a periodic boundary in the vicinal surface in the simulation: if an adatom cross over the upper boundary in the [101: 0] direction, it will reappear in the lower boundary, and it will be the same for a reverse case. In its perpendicular direction, we used a normal periodic boundary condition. The density of steps is a major parameter which is essential to describe the growth processes. The growth mode can be determined from information
CE N IE SC RS CE E A TT RF LE SU
H. Mao et al. / Surface Science 432 (1999) L617–L620
on the density of steps. Especially, if it is compared with the reflection high-energy electron diffraction (RHEED) intensity, some important parameters, such as the energy barrier to hopping, can be determined from the variation of the density of steps. In order to compare our simulation results with some available experimental data [10], we estimated the growth processes on a vicinal surface. In the following simulation shown in Fig. 2, we take the deposition rate to be 1.66 ML per second for the comparison with the experimental results. Two distinct growth modes are noticeable in Fig. 2, depending on the hopping rates. For the hopping rate in the range of 104 to 105 s−1, the density of steps clearly shows an oscillatory behavior, implying that the GaN growth is in a two-dimensional growth mode. However, the oscillatory behaviour in the density of steps with the growth time gradually disappears with the increase in the hopping rate. For the hopping rate D=9×106 s−1; however, there is no oscillation in the density of steps, which implies that the growth mode has been changed to a step flow mode on the vicinal surface. In the RHEED intensity measurements of the GaN growth by MBE, Grandjean et al. found that the growth mode changed from the two-dimensional mode to the step flow mode at a critical temperature T =800°C and at a growth rate of C 1.66ML s−1 [10]. In their experiments the substrate misorientation was less than 2°. With a critical temperature T =800°C and a hopping rate C
Fig. 2. The density of steps on the vicinal surface for different hopping rates are shown. The vicinal surface is misoriented by 2° to the [101: 0] direction from the standard [0001] direction.
L619
D=9×106 s−1, we could roughly estimate the upper limit of the energy barrier to hopping to be 1.43 eV using the Arrhenius expression. There are several factors that affect the precision of the estimated energy barrier to hopping E : the prer factor D , the substrate misorientation and the 0 critical temperature T . The error bar of the energy C barrier to hopping E can be estimated using Eq. r (1). Among the factors mentioned, the prefactor D is the major factor leading to an error. The 0 error due to the prefactor is determined by the substrate temperature and Eq. (2), and it may lead to an uncertainty of 0.2 eV. The uncertainty in the critical temperature is less than 10°C. It only introduces an uncertainty of less than 1%. An uncertainty in the other factors is about the same. Therefore the total uncertainty in the energy barrier to hopping E is less than 0.25 eV. The r upper limit of the energy barrier to hopping is estimated to be 1.43±0.25 eV. This estimated result is very close to an experimental value of 1.45±0.25 eV observed on a Ga-terminated GaN surface by Liu et al. [11]. The activation energy of the surface migration, 1.45±0.25 eV was determined from the recovery behavior of the RHEED intensity[11]. Our model can be further refined. For example, a thermal desorption of a nucleus can be included in the model. If the thermal desorption of a nucleus is considered, the energy barrier to hopping must include two energy terms: a surface contribution E and a contribution from S the nearest neighbor along the surface E . The N values E and E can be regarded as the effective S N energy barriers to hopping. Those values incorporate, in an average way, the fast processes as well as other factors such as the surface reconstruction. According to our previous result, two parameters, E and E /E which will be included in the refined S S N model, are primarily responsible for the temperature of the substrate and the effect of lateral interactions at a given temperature, respectively [9]. It is expected that the values of E and E are S N strongly related to the growth conditions. Therefore our model based upon Monte Carlo simulation and the solid-on-solid model is a useful method to determine the energy barrier to hopping if accurate experimental results are available. In conclusion, we have employed the Monte
SU RF A LE CE TT S ER CIE S N CE
L620
H. Mao et al. / Surface Science 432 (1999) L617–L620
Carlo method and the solid-on-solid model to simulate an epitaxial growth mechanism of the hexagonal structure of GaN semiconductor. The adatoms on the hexagonal structure have more probability to hop freely than those on the cubic structure, resulting in the lower density of islands. The simulation of the growth process on the vicinal surface shows that an epitaxial growth is in a stepflow mode at simulation temperatures above the critical temperature T ; otherwise it is in a twoC dimensional mode. The upper limit of the energy barrier to hopping is estimated to be 1.43±0.25 eV, which is in a good agreement with the reported experimental result.
Acknowledgement This work was supported in part by Korea Science and Engineering Foundation ( KOSEF ).
References [1] S. Nakamura, M. Senoh, N. Iwasa, S. Nagahama, Jpn. J. Appl. Phys. 34 (1995) L797. [2] I. Akasaki, H. Amano, S. Sota, H. Sakai, T. Tanaka, M. Koike, Jpn. J. Appl. Phys. 34 (1995) L1517. [3] S. Nakamura, M. Senoh, S. Nagahama, N. Iwasa, T. Yamada, T. Matsushita, H. Kioku, Y. Sugimoto, Jpn. J. Appl. Phys. 35 (1996) L74. [4] B.W. Lim, Q.C. Chen, J.Y. Yang, M.A. Khan, Appl. Phys. Lett. 68 (1996) 3761. [5] M.A. Khan, Q. Chen, C.J. Sun, J.W. Yang, M.S. Shur, H. Park, Appl. Phys. Lett. 68 (1996) 514. [6 ] J.R. Jenny, R. Kaspi, K.R. Evans, J. Cryst. Growth 175–176 (1997) 89. [7] S. Clarke, D.D. Vvedensky, Phys. Rev. Lett. 58 (1987) 2235. [8] J.D. Weeks, G.H. Gilmer, Adv. Chem. Phys. 40 (1979) 157. [9] H.B. Mao, W. Lu, S.C. Shen, J. Cryst. Growth 151 (1995) 31. [10] N. Gradjean, J. Massies, Appl. Phys. Lett. 71 (1997) 1816. [11] H. Liu, J.G. Kim, M.H. Ludwig, P.M. Park, Appl. Phys. Lett. 71 (1997) 347.