Three-Dimensional Modeling of Basal Plane Dislocations in 4H-SiC

Article pubs.acs.org/crystal

Three-Dimensional Modeling of Basal Plane Dislocations in 4H-SiC Single Crystals Grown by the Physical Vapor Transport Method Bing Gao* and Koichi Kakimoto Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580, Japan ABSTRACT: To effectively reduce basal plane dislocations (BPDs) during SiC physical vapor transport growth, a threedimensional model for tracking the multiplication of BPDs has been developed. The distribution of BPDs inside global crystals has been shown. The effects of the convexity of the growth surface and the cooling rate have been analyzed. The results show that the convexity of the growth surface is unfavorable and can cause a large multiplication of BPDs when the crystal grows. Fast cooling during the cooling process is beneficial for the reduction of BPDs because fast cooling can result in a smaller radial flux at the high-temperature region. In addition, fast cooling can reduce the generation of stacking faults during the cooling process. Therefore, to reduce BPDs and stacking faults, it is better to maintain or reduce the convexity of the growth surface and increase the cooling rate during the cooling process.

1. INTRODUCTION Wide-bandgap silicon carbide (SiC) is a promising semiconductor material for electronic and optoelectronic devices involving high power, high temperature, high frequency, and intense radiation owing to its stable chemical and thermal properties.1,2 Bulk SiC crystals are commonly grown by the physical vapor transport (PVT) method. In bulk SiC crystal growth, significant progress has been achieved in reducing the most damaging defects: micropipes. From 2010 to 2013, researchers at Cree Inc. reduced the micropipe density (MPD) by 90% in 150-mm substrates.3 As the density of micropipes in SiC crystals has been suppressed to a technologically tolerable level, the quality improvement focus has shifted to less severely damaging defects such as dislocations.3 There are several types of dislocations, including basal plane dislocations (BPDs), which are deformation-induced dislocations, and grown-in dislocations (threading edge dislocations (TEDs) and threading screw dislocations (TSDs)). Deformation-induced dislocations usually lie on the primary slip plane, which is the (0001) basal plane. Because BPDs have been associated with an increase in the number of defects observed in epitaxy and as the root cause of VF drift (forward voltage instability) in bipolar devices,4 reduction of the BPD density in SiC crystals has become a major focus in future quality improvement efforts.3 BPDs are mainly generated and multiply during hightemperature processes, such as the crystal growth process, and cooling processes. Thus, optimization of crystal growth and cooling processes could reduce the generation of BPDs. However, study of the optimization process requires a good model that can correctly connect the generation of BPDs to the practical operational conditions and correctly describe the ratedependent plastic deformation process of the SiC crystal. Gao el al.5 extended the Alexander−Haasen model to SiC crystals by including different activation enthalpies in different temperature © 2014 American Chemical Society

regions. Their model showed good agreement between numerical and experimental data for stress−strain curves in a wide temperature range. However, when they applied their model to practical PVT growth of SiC, they did not distinguish between deformation-induced dislocations (BPDs) and grownin dislocations (TEDs and TSDs) and also regarded the von Mises stress as the driving force for dislocation multiplication. It is known that plastic deformation of 4H- and 6H-SiC is mainly caused by BPDs,6−8 and the driving force for BPD multiplication is not the von Mises stress but the resolved shear stress in the primary slip direction, that is, the a/3 < 112̅0 > {0001} direction. Therefore, to correctly track plastic deformation during high-temperature PVT growth, it is necessary to only consider deformation-induced dislocations (i.e., BPDs) and use resolved shear stress along the primary slip directions as the driving force for BPD multiplication. The growth defects, TSD and TED, are neglected during the calculation since they do not affect plastic deformation. In this paper, we improved the model proposed by Gao et al.5 by first resolving thermal stress in the primary slip direction in 4H-SiC and then substituting the resolved shear stress into the three-dimensional (3D) Alexander−Haasen model to obtain the BPDs. The distribution of BPDs inside the crystal and the effect of cooling rate on the final BPD density and stacking faults are also determined.

2. MODELING BASAL PLANE DISLOCATIONS OF 4H-SIC SINGLE CRSYTALS Hexagonal 4H-SiC crystals have three primary slip directions in the (0001) basal plane,9 as shown in Figure 1. The three Received: November 29, 2013 Revised: January 14, 2014 Published: January 17, 2014 1272

dx.doi.org/10.1021/cg401789g | Cryst. Growth Des. 2014, 14, 1272−1278

Crystal Growth & Design

Article

where s is the unit vector of the slip direction in the Cartesian coordinate system, and n is the unit normal vector of the slip plane in the Cartesian coordinate system. After the resolved shear stress is obtained, the basal plane dislocations can be calculated from the following equations:11 dεpl(α)

= Nm(α)v(α)b

(4)

dNm(α) (α)λ (α) (α) = Kτeff Nm v dt

(5)

dt

where the subscript m denotes the mobile dislocation, b is Burgers vector, N is the dislocation density, τeff is the effective stress for dislocation multiplication, K is the multiplication constant, λ is the stress exponential factor, εpl is the creep strain, and ν is the slip velocity of dislocation defined by5,12 ⎛ Q ⎞ (α)p exp⎜ − v(α) = k 0τeff ⎟ ⎝ k bT ⎠

where k0 is a material constant, p is the stress exponential factor, Q is the activation enthalpy, kb is the Boltzmann constant, and T is temperature in Kelvin. The effective stress necessary for dislocation motion is given by12

Figure 1. Hexagonal unit cell. There are two coordinate systems: hexagonal coordinate axes (a1, a2, a3, c) and Cartesian coordinate axes (x, y, z).

(α) (α) τeff = ⟨τresolv − τi(α) − τb(α)⟩

primary slip directions are −a1, −a2, and −a3, which are defined as [2̅110], [12̅10], and [112̅0] in the hexagonal coordinate system and [1̅ 3 0], [1̅ 3 0], and [100] in the Cartesian coordinate system (x, y, z). The resolved shear stress in each slip direction can be calculated by the tensor transformation technique using stress components obtained from a threedimensional (3D) analysis of the stress. However, calculation of the 3D stress requires a large amount of computational time. Thus, a cost-effective method that approximately treats the stress in the axisymmetric case is usually adopted,10 since an ingot has an axisymmetric shape and the thermal condition during the ingot cooling process is also axisymmetric. The radial stress σrr, axial stress σzz, azimuthal stress σφφ, and shear stress σrz can be obtained from an axisymmetric thermal stress analysis and then converted into the resolved shear stresses τ(α) for the three primary slip systems using the tensor transformation technique. The basic procedure is the following. First, change the stress tensor from the cylindrical coordinate system (r, φ, z) to the Cartesian coordinate system (x, y, z), where z is the [0001] direction:

σcart

⎡ σrr 0 σrz ⎤ ⎡ σxx σxy σxz ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ σxy σyy σyz ⎥ = A⎢ 0 σφφ 0 ⎥AT ⎢ ⎥ ⎢σ σ σ ⎥ ⎣ xz yz zz ⎦ ⎣ σrz 0 σzz ⎦

τi(α) = μb

∑ aαβNm(β) (8)

β

τb(α) = μb ∑ Aαβ Nm(β) (9)

β

Owing to the equivalence of the three slip directions, the values of Aαβ and ααβ are taken as fixed values: ααβ = 0.09, and Aαβ|α=β = 0.125 or Aαβ|α≠β = 0.0625. After the dislocation densities and creep strains are calculated for all of the slip directions, the total dislocation density Nm and total creep strain tensor εpl in the Cartesian coordinate system can be expressed as 3

Nm =

∑ Nm(α) α=1

(1)

3

εpl =

∑ εpl(α) α=1

(10)

1 (α) (α) (n ⊗ m(α) + m(α) ⊗ n(α))sign(τresolv ) 2 (11)

Table 1 gives all the parameters for the calculations.

3. NUMERICAL RESULTS For the configurations and size of the simulated PVT furnace, refer to refs 5 and 13. The growth chamber has a shape of frustum, which has a top diameter of 50 mm, a bottom diameter of 100 mm, and a height of 43 mm. For the numerical simulations of the crystal growth and cooling process, the

(2)

The resolved shear stress along the α slip direction can be easily obtained by (α) τresolv = s(α)·σcart·n(α)

(7)

where τ(α) is the stress required to overcome short-range i obstacles, and τ(α) is the internal long-range elastic stress b generated by mobile dislocations.12 The angle brackets ⟨ ⟩ denote a non-negative value; that is, ⟨x⟩ = x if x > 0, and 0 if x ≤ 0. The short- and long-range interactions are given by12

where φ is the angle between the [100] direction and any vector in the basal plane and ⎡ cosφ −sinφ 0 ⎤ ⎢ ⎥ A = ⎢ sinφ cosφ 0 ⎥ ⎢⎣ 0 0 1 ⎥⎦

(6)

(3) 1273



Article

Table 1. Parameters for the Dislocation Calculations symbol

description

b K k0 Q

Burger’s vector multiplication constant material constant5 activation enthalpy5

λ

stress exponential factor5

p kb μ Nm τeff ν s n εpl εpl α,β T

stress exponential factor5 Boltzmann constant Shear modulus mobile dislocation density effective stress slip velocity of dislocation unit vector of slip direction unit normal vector of slip plane plastic strain in slip direction plastic strain tensor in Cartesian system slip direction temperature

value and/or units 3.073 × 10−10 (m) 7.0 × 10−5 8.5 × 10−15 3.3 eV at T > 1000 °C 2.6 eV at T < 1000 °C 1.1 at T > 1000 °C 0.6 at T < 1000 °C 2.8 8.615 × 10−5 (eV/K) 7.992 × 1010 (Pa) m−2 Pa m/s / / / / / K

Figure 2. 3D view of the BPD distribution at room temperature for 10 h cooling time.

104 cm−2. Thus, the present crystal exhibits a low dislocation density in most of the crystal except for the convex part near the bottom of the crystal. From the top to the bottom of crystal, the BPD density gradually decreases and then gradually increases. The high BPD density at the top of crystal is understandable since there is a large cooling flux exiting from the top of crystal. The high BPD density at the bottom of crystal is due to a convex interface as explained in the following section. In the middle part, the curvature of interface is quite flat, and therefore, a low BPD density exists there. 3.2. Analysis of the Largest BPD Generation. As shown in Figure 3a2−d2, there are large BPDs generated in the convex part of the bottom of the crystal. Since large BPDs are not desirable, it is essential to analyze the reason for the generation of large BPD and avoid it during the crystal growth process. BPDs can be generated during both the growth and cooling processes. It is important to determine which process is dominant for the largest BPD generation at the convex crystal bottom. Figure 4 shows the BPD distribution at the end of growth process, that is, before the cooling process. It shows that most of the BPDs are already generated during the growth process, although Figure 4 shows a slightly smaller area of the BPD distribution at the crystal bottom and less BPD generation at the crystal side compared with Figure 2. Therefore, it is the growth process that caused the large BPD generation at the bottom of the crystal. To determine the main reason for the generation of large BPDs, it is necessary to determine the time that the BPDs are generated during the growth process. Figure 5 shows the BPD distribution at growth times of 35, 42.5, 50, and 57.5 h. It shows that the large BPDs were mainly generated between 42.5 h (Figure 5b) and 50 h (Figure 5c). After the generation time was obtained, the temperature distributions for these growth times are shown in Figure 6. The arrows in Figure 6 show the location of the maximum BPDs. From Figure 6a−d, the temperature at the bottom of the crystal gradually increases and results in a large dislocation mobile velocity. Furthermore, from Figure 6a−d, the vertical distances from the center to the edge of the bottom of crystal are 3.6, 4.1, 4.4, and 4.6 mm, respectively, and thus the convexity gradually increases and causes a large thermal stress component. For on-axis or very small off-axis growth of SiC crystals, the dominant stress component for BPD multiplication is σrz.9 Figure 7 shows the distribution of σrz for the temperature in Figure 6. It shows that there is a large value of σrz along the bottom of the crystal. Therefore, to control the BPD density to a constant small value

power of the heating system was adjusted by controlling the intensity of the electric current inside the coils. The square of the electrical current intensity was first kept constant during the growth time (60 h) and then linearly decreased during two cooling times (1 and 10 h). The axisymmetric temperature distribution inside the crystal is shown in ref 5. After the temperature field was obtained, the coupled calculations of stress and dislocation were performed. In the stress calculations, free boundary conditions are applied on the side boundary because our furnace configurations can effectively avoid poly-SiC formed surrounding the single crystal and avoid contact between crystal and crucible. The boundary conditions between the seed and the lid of the crucible are complicated and can be taken as three kinds of boundary conditions as explained in the paper.14 In the present case, only free boundary conditions are applied because of the freestanding setup of the seed in our experiments. The calculation results were given in four parts: the first part showed the distribution of the BPD density inside the global crystal; the second part gave an analysis for the largest BPD generation; and the third and fourth parts showed the effects of the cooling process on the generation of BPDs and stacking faults, respectively. 3.1. Distribution of BPD Density at Room Temperature. The 3D BPD density at room temperature for the 10 h cooling time is shown in Figure 2. The crystal grows from the top to the bottom. The surface of the bottom of the crystal is slightly convex. The order of the BPD density at the bottom of the crystal is 104 cm−2, and it shows 6-fold symmetry on the basal plane, which is consistent with the experimental data.15 The distribution of the BPD density inside the global crystal is shown by eight slices from the top to the bottom in Figure 3. In Figure 3a1−d1, the distance between slices is 4.7 mm, and the slice d1 is located at the same height as the edge of the convex surface of the bottom. The slices from a2−d2 all belong to the convex part near the bottom of the crystal, and the distance between slices is 1.1 mm. It can be seen that most of the crystal has low BPD density (1750 K). Therefore, the reason 1276



Article

Figure 9. (a) Variation of the maximum radial flux for 1 and 10 h cooling times, and (b) variation of the maximum axial flux gradient for 1 and 10 h cooling times.

Figure 10. Increase in the BPD density from 1000 °C to room temperature during the cooling process: (a) 1 h cooling time and (b) 10 h cooling time.

which the mobility of BPDs becomes low. Thus, the increase of the radial cooling flux does not contribute so much to the generation of BPDs. Figure 9 also shows that the generation of dislocations mainly occurs in the high-temperature region. Thus, good control in the high-temperature region is crucial to reducing the dislocations. 3.4. Effect of Cooling Rate on Stacking Faults. The deformation-induced BPDs in hexagonal SiC can be easily dissociated into widely separated partial dislocations. The leading partial (Si core) moves more easily than the trailing partial (C core). When the temperature is below 1000 °C, it is possible that the trailing partial may not move at all, and the leading partial will glide along, leaving behind a stacking fault.18−20 When the temperature is higher than 1000 °C, and the stacking faults contributed by the preferential movement of Si partial dislocations are small, because both Si-core partial and C-core partial can move at high temperature, and the movement of trailing partial may eliminate the stacking faults left by the leading partial. Therefore, when the temperature is below 1000 °C during the cooling process, the region of the increase in the dislocations might correspond to the region of stacking faults. To determine the possible region of stacking faults, the increase of dislocations was obtained using the final dislocations at room temperature minus the dislocations at 1000 °C for every point. Figure 10 outlines the increase in the dislocations

1 h cooling results in less BPDs must be because of the smaller radial flux in the high-temperature region. Why does fast cooling (1 h) cause less radial flux in the hightemperature region? When the crystal growth process finishes, the crucible temperature is much higher than the crystal temperature, which causes heating flux. When the cooling process begins, fast cooling causes a rapid decrease of the crucible temperature and thus causes a rapid decrease of heating flux. This is the main reason that fast cooling results in a rapid decrease of radial flux in the high-temperature region, which corresponds to stage 1 for the 1 h case in Figure 9a. When the crucible temperature is equal to the crystal temperature, the radial flux becomes zero. (Zero radial flux does not appear in Figure 9a because Figure 9a only shows the maximum value of radial flux among all of points inside crystal.) However, when the crucible temperature is lower than the crystal temperature, radial heating flux becomes the cooling flux. The rapid decrease of the crucible temperature causes a rapidly increasing cooling flux, which corresponds to stage 2 for the 1 h case in Figure 9a. After cooling for some time, the total heat inside crystal is reduced, and the cooling flux becomes weak and gradually decreases, which corresponds to stage 3 for the 1 h case in Figure 9a. Therefore, fast cooling is beneficial to reduce the generation of BPDs if the heat inside the crystal can be easily dissipated through the crucible, like in the PVT case. Figure 9a shows that when the radial cooling flux starts to increase, the temperature is already decreased to 1910 K, at 1277



Article

from 1000 °C to room temperature in a crystal for cooling times of 1 and 10 h. It shows that a long cooling time results in an increase in BPD density, which means that a long cooling time results in more stacking faults. Therefore, to reduce the stacking faults, it is better to use fast cooling.

(15) Sakwe, S. A.; Stockmeier, M.; Hens, P.; Müller, R.; Queren, D.; Kunecke, U.; Konias, K.; Hock, R.; Magerl, A.; Pons, M.; Winnacker, A.; Wellmann, P. Silicon Carbide, Vol. 1: Growth, Defects, and Novel Applications; Friedrichs, P., Kimoto, T., Ley, L., Pensl, G., Eds.; WileyVCH: New York, 2011; pp 1−31. (16) Gao, B.; Nakano, S.; Harada, H.; Miyamura, Y.; Kakimoto, K. Cryst. Growth Des. 2013, 13, 2661−2669. (17) Gao, B.; Kakimoto, K. J. Cryst. Growth 2013, 384, 13−20. (18) Lara, A.; Muñoz, A.; Castillo-Rodríguez, M.; DomínguezRodríguez, A. Ceram. Int. 2012, 38, 1381−1390. (19) Pirouz, P.; Demenet, J. L.; Hong, M. H. Phil. Mag. A 2001, 81, 1207−1227. (20) Samant, A. V.; Hong, M. H.; Pirouz, P. Phys. Status Solidi B 2000, 222, 75−93.

4. CONCLUSIONS Three-dimensional modeling of basal plane dislocations in 4HSiC single crystal grown by the PVT method has been performed. The results show that the BPD density has 6-fold symmetry inside the crystal. Although the large convexity of the growth surface is beneficial for reducing grown-in dislocations, it causes large basal plane dislocations when the crystal grows up. Thus, during the global crystal growth process, it is better to maintain or reduce the convexity of growth. The cooling process has some effect on the generation of BPDs. To reduce the generation of BPDs, it is better to use fast cooling, because fast cooling can result in lower radial flux in the hightemperature region. In addition, fast cooling can generate smaller stacking faults during the cooling process compared with slow cooling. Therefore, to reduce BPDs and stacking faults, it is better to reduce the convexity of the growth surface and increase the cooling rate during the cooling process.

■

AUTHOR INFORMATION

Corresponding Author

*Tel: +81-92-583-7744. Fax: +81-92-583-7743. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was partly supported by the Japan Society for the Promotion of Science (Grant No. 24360012). REFERENCES

(1) Wang, X. L.; Cai, D.; Zhang, H. Int. J. Heat Mass Transfer 2007, 50, 1221−1230. (2) Wang, X. L.; Cai, D.; Zhang, H. J. Cryst. Growth 2007, 305, 122− 132. (3) Powell, A. R.; Leonard, R. L.; Khlebnikov, Y.; Deyneka, E.; Mckay, M.; Sumakeris, J. J.; Tsvetkov, V.; Balkas, E. Abstracts of International Conference on Silicon Carbide and Related Materials 2013; September 29−October 4, 2013, Miyazaki, Japan; The Japan Society of Applied Physics: Tokyo, Japan, 2013; No. 91. (4) Lendenman, H.; Dahlquist, F.; Johansson, N.; Soderholm, R.; Nilsson, P. A.; Bergman, J. P.; Skytt, P. Mater. Sci. Forum 2001, 353− 356, 727−730. (5) Gao, B.; Kakimoto, K. J. Cryst. Growth 2014, 386, 215−219. (6) Fujita, S.; Maeda, K.; Hyodo, S. Philos. Mag. A 1987, 55, 203− 215. (7) Suematsu, H.; Suzuki, T.; Iseki, T.; Mori, T. J. Am. Ceram. Soc. 1991, 74, 173−178. (8) Zhang, M.; Hobgood, H. M.; Pirouz, P. Mater. Sci. Forum 2004, 457−460, 371−374. (9) Bottcher, K.; Cliffe, K. A. J. Cryst. Growth 2007, 303, 310−313. (10) Miyazaki, N.; Kuroda, Y.; Sakaguchi, M. J. Cryst. Growth 2000, 218, 221−231. (11) Orowan, E. Philos. Trans. R. Soc. London Ser. A 1940, 52, 8−22. (12) Cochard, J.; Yonenaga, I.; Gouttebroze, S.; M’Hamdi, M.; Zhang, Z. L. J. Appl. Phys. 2010, 107, 033512. (13) Gao, B.; Chen, X. J.; Nakano, S.; Nishizawa, S.; Kakimoto, K. J. Cryst. Growth 2010, 312, 3349−3355. (14) Ma, R. H.; Zhang, H.; Prasad, V.; Dudley, M. Cryst. Growth Des. 2002, 2, 213−220. 1278


Three-Dimensional Modeling of Basal Plane Dislocations in 4H-SiC

Recommend Documents