Journal of Crystal Growth 351 (2012) 32–36
Contents lists available at SciVerse ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Molecular dynamics simulation of diffusion behavior of N atoms on the growth surface in GaN solution growth Takahiro Kawamura a,n, Yoshihiro Kangawa b, Koichi Kakimoto b, Yasuyuki Suzuki a a b
Graduate School of Engineering, Mie University, 1577 Kurimamachiya-cho, Tsu, Mie 514-8507, Japan Research Institute of Applied Physics, Kyushu University, 6-1 Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan
a r t i c l e i n f o
a b s t r a c t
Article history: Received 24 October 2011 Received in revised form 21 February 2012 Accepted 14 April 2012 Communicated by M. Uwaha Available online 23 April 2012
In this study, we simulated the solution growth of gallium nitride (GaN) and investigated the diffusion behavior of nitrogen (N) atoms on growth surfaces by molecular dynamics simulation. The simulation showed that the Ga-face grew flatter than the N-face. Comparing the diffusion coefficients of a N atom on Ga- and N-faces, the values on the Ga-face were about 3.5 times larger than those on the N-face. & 2012 Elsevier B.V. All rights reserved.
Keywords: A1. Computer simulation A1. Diffusion A1. Molecular dynamics A2. Growth from solutions B1. GaN
1. Introduction Gallium nitride (GaN)-based devices are expected to comprise the next generation of light-emitting, high-power and energysaving devices because of their physical properties, such as wide bandgap, high breakdown field, high saturated electron-drift velocity [1] and high thermal conductivity [2,3]. These devices are generally fabricated on sapphire or silicon carbide (SiC) substrates, but because of lattice mismatch and differences in thermal expansion coefficients, epitaxial layers contain lattice defects that result in deteriorated device performance. The most promising approach to solving the problem is fabrication on a GaN substrate. However, because it is difficult to grow highquality, large bulk GaN single crystals, the GaN substrate is expensive and cannot be used easily. Therefore, many researchers have attempted to develop growth techniques for the production of large bulk GaN single crystals. Bulk GaN single crystals have been grown by several methods, such as hydride vapor phase epitaxy (HVPE) [4,5], ammonothermal growth [6,7], high-pressure solution growth (HPSG) [8,9] and Na-flux solution growth [10,11]. The HVPE method results in a faster growth rate (above 100 mm=h) than others, but its high production cost and by-product GaCl4 are issues to be resolved.
n
Corresponding author. Tel.: þ81 59 231 9364. E-mail address:
[email protected] (T. Kawamura).
0022-0248/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jcrysgro.2012.04.022
The HPSG method can grow high-quality crystal, but it is difficult to grow large crystals. The ammonothermal and Na-flux solution growth methods succeeded in producing large bulk crystals [6,11], but the increase in growth rate and enlargement of the crystal size are common concerns about these liquid phase epitaxy methods at present. The further improvement of crystal growth techniques requires a detailed understanding of the crystal growth process on the atomic scale. Molecular dynamics (MD) simulation is a powerful approach for investigating the crystal growth mechanism, because it gives us atomic-scale information. Many studies on crystal growth simulation for semiconductors using the MD method have been reported [12–16]. We previously simulated solution growth of GaN using the MD method to investigate the difference in growth mechanisms on Ga- and N-faces [18]. The simulation result showed that the growth surface on the Ga-face grew more flatly than that on the N-face. We considered the difference in the surface morphology on the Ga- and N-faces to be due to the diffusion behavior of nitrogen (N) atoms on the growth surface. It was because that the N concentration in Ga solution is extremely low, the diffusion behavior of N atoms appears to strongly affect its growth process and crystal quality. However, we were not able to show any objective evidence of this conclusion at that time. Therefore, in the present study we examined the diffusion behavior of N atoms on the growth surfaces and investigated the influence on surface morphology. First we show typical simulation results of the
T. Kawamura et al. / Journal of Crystal Growth 351 (2012) 32–36
solution growth of GaN on Ga- and N-faces. Second, we calculate the diffusion coefficients of a N atom on the growth surface of Gaand N-faces.
2. Simulation methods 2.1. Molecular dynamics We employed a classical MD method based on the law of Newton’s equation of motion. The Brenner-type potential function was used in our simulation [19]. The potential functions are given by XX F¼ ½V R ðr ij ÞBnij V A ðr ij Þ, ð1Þ i
ioj
V R ðr ij Þ ¼ f ðr ij Þ
qffiffiffiffiffiffiffiffi Dij expð 2Sij bij ðr ij Rij ÞÞ, Sij 1
sffiffiffiffiffi ! Dij Sij 2 exp b ðr R Þ : V A ðr ij Þ ¼ f ðr ij Þ Sij ij ij ij Sij 1
ð2Þ
ð3Þ
Eqs. (2) and (3) are repulsive and attractive interactions, respectively. The Bij is a many-body effect that represents the contribution from atom k around atom i or j. The parameters for GaN were fitted by Nord et al. [20], and they showed that the potential gives a reasonable description of the melting behavior and solubility of N.
Fig. 1 shows a simulation cell of solution growth, which consists of crystal and solution parts. The size of the simulation cell is 1.9 1.6 5.0 nm3. The crystal part is composed of 216 Ga and 216 N atoms. The solution part is composed of 476 Ga and up to 100 N atoms. We inserted N atoms one by one at the center of the solution region every 60 ps during the growth simulation. Under the simulation condition, the N concentration was extremely higher than the practical case [17]. However, it is impracticable to simulate solution growth under the realistic N concentration condition because it needs very large simulation cell. We decided the number of atoms in the simulation cell in accordance with our computing environment. The simulation procedure of the solution growth is described in Ref. [18] in detail. The simulation conditions are as follows. We used the Gear algorithm for numerical integration. The time step was 0.3 fs, and the total simulation time for solution growth and diffusion behavior was 15 and 3 ns, respectively. The temperature was controlled by the Nose–Hoover method [21,22]. The periodic boundary condition was used in the direction perpendicular to the crystal/solution interface. The specular reflection boundary condition was used at the top surface of the simulation cell. The atoms in the bottom layer of the crystal part were fixed. 2.2. Diffusion coefficient We used a similar simulation cell to that shown in Fig. 1 for the simulation of diffusion behavior. We initially set a N atom in the interface between the crystal and solution parts. In this simulation, we intended to simulate diffusion behavior of N atoms on the crystal surface (crystal/solution interface). Thus we applied the specular reflection boundary condition to only the N atom to restrict the region that the N atom migrates. Fig. 2 shows the concept of the restriction. The N atom can migrate in the region from the crystal surface to the height that the specular reflection boundary is set. We estimated the influence of the height on diffusion behavior. The result is shown in the last of Section 3.2. We examined the diffusion trajectory of the N atom and calculated diffusion coefficient D using the following equation: D¼
Fig. 1. Simulation cell.
33
1 2 /9rðt 0 þ tÞrðt 0 Þ9 S, 6t
ð4Þ
where vector r is the atomic coordinate, t is the simulation time, 2 r(t0) is the initial coordinate and 9rðt 0 þ tÞrðt 0 Þ9 is the mean square displacement (MSD). Fig. 3 shows the result of MSD of the N atom on the Ga-face simulated at 2500 K. The diffusion coefficient is estimated from the slope of the MSD. However, when the migration frequency is small and the values of the MSD discontinuously changes like the result, it is difficult to estimate the slope. In this study, because the number of target N atom was
Fig. 2. Restriction of the N migration on the crystal surface.
34
T. Kawamura et al. / Journal of Crystal Growth 351 (2012) 32–36
Fig. 5. Averaged mean square displacement.
Fig. 3. Normal mean square displacement.
Fig. 4. Schematic diagram of sampling data of mean square displacement.
one, the shortage of the sampling number also affected the result. Therefore, we recorded the data of MSD many times changing the initial position r(t0) during the simulation and finally averaged them to obtain an averaged result of MSD that the slope was estimative. Fig. 4 shows a schematic diagram of sampling data of MSD. In this case, an averaged diffusion coefficient is calculated by D¼
n X m 1X 1 9rðt 0 ði1ÞDt þ jDtÞ n i ¼ 1 j ¼ 1 6t s 2
rðt 0 ði1ÞDtÞ9 ,
ð5Þ
ð6Þ
where Dt is step interval recording the data of r, ts ¼mDt and n is the number of averaging. The values of each parameter were Dt ¼ 1 103 steps, ts ¼1 106 steps and n ¼9 103. Fig. 5 shows the averaged MSD. The results shown in Figs. 3 and 5 are obtained from the same simulation. Because the averaged MSD shows smoother result than the normal MSD (Fig. 3), we can easily estimate the slope of MSD from the averaged MSD. The activation energy for migration E was estimated from temperature dependence of diffusion coefficient by E D ¼ D0 exp , ð7Þ kB T where D0 and kB are the pre-exponential factor and the Boltzman constant, respectively.
Fig. 6. Snapshots of GaN solution growth on the Ga-face.
3. Results and discussion 3.1. Solution growth simulation on Ga- and N-faces Figs. 6 and 7 show snapshots of GaN solution growth on Ga- and N-faces, respectively. The temperature was set to 2500 K, which is close to the melting point of GaN [23], to prevent crystallization in solution. However, because the N concentration was extremely high under the simulation condition, spontaneous nucleation occurred in the solution. After the nucleation, crystal growth on the substrate surface did not develop because most of the inserted N atoms adsorbed on the nucleation. We observed the growth processes on the Ga- and N-faces until the crystal growth had stopped. We can see in Fig. 6 that the growth surface on the Ga-face was flat. We observed that N atoms those reached on the growth surface diffused in the solid–liquid interface. Because the N atoms unstably bonded to the crystal surface by one bond, the N atoms appeared to easily diffuse in the crystal/solution interface. On the other hand, it is found in Fig. 7 that the N-face was rough. The snapshot after 4.5 ns shows that only part of the crystal surface
T. Kawamura et al. / Journal of Crystal Growth 351 (2012) 32–36
35
Fig. 8. Temperature dependence of diffusion coefficients of a N atom on the Gaand N-faces. Fig. 7. Snapshots of GaN solution growth on the N-face.
developed. After 6 ns, a GaN cluster that generated in the solution bonded to the prominent part and the N-face became more rough. The result that the Ga-face grew flatter than the N-face was consistent with the experimental results by the HPSG method [8,9]. We speculated that the difference in the surface morphology on the Ga- and N-faces was due to the different diffusion behaviors of N atoms on the growth surfaces. Because the surface diffusion length of N atom on the Ga-face might be longer than that on the N-face, the Ga-face grew flatter than the N-face. In order to confirm this speculation, we calculated diffusion coefficient of a N atom in the growth surface (crystal/solution interface) and compared the values on the Ga- and N-faces, as described in the next subsection.
Table 1 Activation energy for migration E and pre-exponential factor D0 calculated from the results shown in Fig. 8.
Ga-face N-face Ga(bulk) N(bulk)
E (eV)
D0 (m2/s)
1.59 7 0.18 2.27 7 0.26 0.37 70.02 0.47 70.12
4.02 10 7 2.50 10 6 7.67 10 8 3.78 10 8
3.2. Diffusion coefficients of a N atom on Ga- and N-faces Fig. 8 shows the Arrhenius plot of temperature dependence of the diffusion coefficients of a N atom in the temperature range from 2000 K to 3000 K. Diffusion coefficients of Ga and N atoms in bulk Ga solution are also shown for reference. The bulk solution simulations were carried out using the simulation cell consisted of 1000 Ga and 5 N atoms. We did not find an experimentally measured diffusion coefficient of N atoms in Ga solution in the literature, nor did we find a Ga self-diffusion coefficient in the same temperature range as our simulation. The calculated value of the self-diffusion coefficient of Ga was reported to be 3.2 10 8 m2/s at 2500 K [24], which was comparable with the value of 1.4 10 8 m2/s in our simulation, confirming the validity of our results to some extent. It is found in Fig. 8 that the values on the Ga-face were higher than those on the N-face in the temperature range. The diffusion coefficient of the N atom on the Ga-face was about 3.5 times larger on average than that on the N-face. The activation energies and pre-exponential factors for the Ga- and N-faces and the Ga and N atoms in the bulk solution were obtained from the results shown in Fig. 8 and they are indicated in Table 1. The value of activation energy for the Ga-face was smaller than that for N-face. These results show that N atoms diffuse more easily on Ga-face than that on N-face. The result agrees well with the speculation obtained from the results of solution growth simulation. Finally, we investigated the influence of the height of the region that a N atom migrates on the diffusion behavior. This
Fig. 9. Dependence of diffusion coefficient of a N atom on the height of the region that the N atom migrates.
investigation aimed to check whether the simulation method that was employed for calculating diffusion coefficients in the crystal/ solution interface affected the diffusion behavior or not. Fig. 9 shows the dependence of diffusion coefficient of a N atom on height of the diffusion region. We simulated under Ga- and N-face conditions. It was found that the values of diffusion coefficient increased with increasing the height of the diffusion region. It was because that when the height of the diffusion region was high, the N atom sometimes migrated in the solution and it is shown in Fig. 8 that the diffusion coefficients of the N atoms in the solution were larger than those on the crystal surface. However, the values on Ga-face were always larger than those on N-face in the height
36
T. Kawamura et al. / Journal of Crystal Growth 351 (2012) 32–36
˚ Therefore, we concluded that the simulation range from 3 to 8 A. method did not affect the results of the diffusion behavior on the crystal surface.
4. Conclusions We simulated the solution growth of GaN and investigated the diffusion behavior of N atoms on the growth surfaces by the MD simulation. The Ga-face grew more flatly than the N-face. The values of diffusion coefficients of N atoms on the Ga-face were about 3.5 times larger than those on the N-face. Therefore, we conclude that because the surface diffusion length of N atoms on the Ga-face was longer than that on the N-face, the Ga-face grew more flatly than the N-face.
Acknowledgments This work was supported in part by the Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu University. References [1] W.J. Choyke, H. Matsunami, G. Pensl, Silicon Carbide: Recent Major Advances, Springer, Berlin, 2004, pp. 737. [2] W. Liu, A.A. Balandin, C. Lee, H.-Y. Lee, Physica Status Solidi A 202 (2005) R135. [3] H. Shibata, Y. Waseda, H. Ohta, K. Kiyomi, K. Shimoyama, K. Fujito, H. Nagaoka, Y. Kagamitani, R. Simura, T. Fukuda, Materials Transactions 48 (2007) 2782.
[4] K. Fujio, S. Kubo, H. Nagaoka, T. Mochizuki, H. Namita, S. Nagao, Journal of Crystal Growth 311 (2009) 3011. [5] E. Richter, U. Zeimer, S. Hagedorn, M. Wagner, F. Brunner, M. Weyers, ¨ G. Trankle, Journal of Crystal Growth 312 (2010) 2537. [6] T. Hashimoto, F. Wu, J.S. Speck, S. Nakamura, Journal of Crystal Growth 310 (2008) 3907. [7] R. Dwilin´ski, R. Doradzin´ski, J. Garczyn´ski, L. Sierzputowski, R. Kucharski, M. Zajac, M. Rudzin´ski, R. Kudrawiec, W. Strupin´ski, Physica Status Solidi A 208 (2011) 1489. [8] M. Boc´kowski, I. Grzegory, S. Krukowski, B. Łucznik, Z. Romanowski, M. Wro´blewski, J. Borysiuk, J. Weyher, P. Hageman, S. Porowski, Journal of Crystal Growth 246 (2002) 194. [9] M. Boc´kowski, P. Strak, I. Grzegory, B. Łucznik, S. Porowski, Journal of Crystal Growth 310 (2008) 3924. [10] T. Yamada, H. Yamane, Y. Yao, M. Yokoyama, T. Sekiguchi, Materials Research Bulletin 44 (2009) 594. [11] M. Imade, Y. Hirabayashi, Y. Konishi, H. Ukegawa, N. Miyoshi, M. Yoshimura, T. Sasaki, Y. Kitaoka, Y. Mori, Applied Physics Express 3 (2010) 075501. [12] T. Motooka, K. Nisihira, S. Munetoh, K. Moriguchi, A. Shintani, Physical Review B 61 (2000) 8537. [13] X.W. Zhou, D.A. Murdick, B. Gillespie, H.N.G. Wadley, Physical Review B 73 (2006) 045337. [14] F. Gao, Y. Zhang, R. Devanathan, M. Posselt, W.J. Weber, Nuclear Instrumentation and Methods in Physics Research B 255 (2007) 136. [15] C. Tang, L. Meng, H. Xiao, J. Zhong, Journal of Applied Physics 103 (2008) 063505. [16] B.A. Gillespie, H.N.G. Wadley, Journal of Crystal Growth 311 (2009) 3195. [17] M. Morishita, F. Kawamura, M. Kawahara, M. Yoshimura, Y. Mori, T. Sasaki, Journal of Crystal Growth 284 (2005) 91. [18] T. Kawamura, Y. Kangawa, K. Kakimoto, S. Kotake, Y. Suzuki, Japanese Journal of Applied Physics 51 (2012) 01AF06. [19] D.W. Brenner, Physical Review B 42 (1990) 9458. [20] J. Nord, K. Albe, P. Erhart, K. Nordlund, Journal of Physics: Condensed Matter 15 (2003) 5649. [21] S. Nose, Journal of Chemical Physics 81 (1984) 511. [22] W.G. Hoover, Physical Review A 31 (1985) 1695. [23] J.A. Van Vechten, Physical Review B 7 (1973) 1479. [24] D.K. Belashchenko, O.I. Ostrovskii, Russian Journal of Physical Chemistry 80 (2006) 509.