- Email: [email protected]
Vacancy behavior in Czochralski silicon growth
ARTICLE IN PRESS Journal of Crystal Growth 311 (2009) 3592–3597
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Vacancy behavior in Czochralski silicon growth Sang Hun Lee a,b,, Jeong Won Kang c, Young Ho Hong a,d, Hyun Jung Oh a, Do Hyun Kim b a
Simulation Part, LG Siltron, Imsoodong, Gumi, Gyeongbuk 730-724, Republic of Korea Department of Chemical and Biomolecular Engineering, KAIST, Daejeon 305-701, Republic of Korea Department of Computer Engineering, Chungju National University, Chungju 380-702, Republic of Korea d Department of Ceramic Engineering Hanyang University, Seoul 133-791, Republic of Korea b c
a r t i c l e in f o
a b s t r a c t
Article history: Received 3 October 2008 Received in revised form 28 April 2009 Accepted 29 April 2009 Communicated by K.W. Benz Available online 14 May 2009
The vacancy defect behavior in silicon crystal growth has been investigated by kinetic Monte Carlo and continuum simulations. The vacancy concentration distributions in silicon crystal were obtained from continuum model simulations corresponding to experimental conditions and then the vacancy clusters distribution obtained from kinetic lattice Monte Carlo simulations. At the temperature above 1470 K, the diffusivity was almost constant because the clusters were not formed. In the clustering temperature region (1370–1270 K), the larger clusters were generated at the higher temperature, the smaller clusters were generated at the lower temperature. While the vacancy concentration led to increase in the number of clusters, the mean size of clusters was irrespective of the vacancy concentrations. Clustering phenomena were very susceptible to temperature and vacancy concentrations. The total number of vacancy clusters was linearly proportion to the crystal pull rate. The radial distribution of clusters obtained from multi-scale simulations was in good agreement with the distribution of voids in the experimental data. & 2009 Elsevier B.V. All rights reserved.
PACS: 61.50.Ah 61.72.Cc 61.72.jd Keywords: A1. Computer simulation A1. Diffusion A1. Kinetic lattice Monte Carlo A1. Point defects B2. Semiconducting silicon
1. Introduction The formation of defects in a silicon crystal during the growth process strongly affects the gate oxide integrity (GOI) characteristics of metal-oxide-semiconductor (MOS) devices [1]. Grown-in defects are formed by the aggregation of excess vacancies or selfinterstitials, and are thought to be related to the distribution of the point defect concentration during the crystal growth process. There are many theoretical and experimental studies for defect formation. Voids are formed by vacancy aggregation [2]. The type of grown-in defect depends on the growth rate. Voronkov [3] showed that the type of grown-in defect depends on the ratio V/G (where V is growth rate and G is axial temperature gradient of the crystal) and that there is a critical ratio that defines the defect type. Such studies prompted comprehensive modeling of point defect dynamics, and many researches have been developing and testing all scale simulation models including quantum mechanical calculations. Habu et al. [4,5] proposed models to describe the
Corresponding author at: Simulation Part, LG Siltron, Imsoodong, Gumi, Gyeongbuk 730-724, Republic of Korea. Tel.: +82 54 713 3799; fax: +82 54 713 3518. E-mail address: [email protected] (S.H. Lee).
0022-0248/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2009.04.044
dominant point defect species, and their model realistically describes the two-dimensional distribution of point defects in Czochralski (CZ) crystals. Ebe [6] performed analytical calculations on the behavior of point defects in the growing silicon crystal. Several theoretical descriptions of the vacancy cluster formation process in silicon exist. Bongiorno and Colombo [7] used the Stillinger–Weber potential [8] and molecular dynamics simulations to demonstrate that strain fields due to the distortion in the silicon lattice near vacancies affect the capture radius in different directions from a vacancy cluster. La Magna et al. [9] proposed a kinetic lattice Monte Carlo (KLMC) model for the description of vacancy diffusion and clustering phenomena. Colombo [10] described native defects and their interaction in silicon using several KLMC models. In this work, the interaction of the point defects in silicon crystal growth has been investigated by the multi-scale method combining KLMC and continuum models. Such a strategy is based on a coherent development of connected simulation tools with different accessibility domains. The continuum model can simulate whole processes but cannot consider atomic interactions directly whereas the KLMC model can simulate dynamic atomic interactions but the analyzing area is limited to nanometer scales. The multi-scale method is expected to play an important role in
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
the investigation of micro defects formation and distribution. In this work, we present the simulation results together with the experiment results corresponding to each simulation to validate our simulations.
2. Simulation methods Fig. 1 shows the process of multi-scale method, which consists of two types of simulation module for analyzing void generation with variable pull rate. Since the continuum model is useful for analyzing the whole system, point defect distributions of the whole crystal are obtained from continuum model simulation. However, this could not describe the aggregate of vacancies [11,12]. The KLMC model makes up for the weak point of the continuum model. This model is based on interactions between individual defects. These compliments are a good approach to understand void generation. Multi-scale method is a strategy to connect different scale simulation modules. We tried to analyze the real condition of industrial processes like changing pull rate. The continuum models for point defect dynamics have been developed by considering diffusion and recombination of defect species like vacancies, self-interstitials and impurities such as oxygen. We consider only native point defects in this continuum model. The concentrations of point defects during the Czochralski process are obtained by solving Eqs. (1) and (2) as follows: DV;I rC V;I Q V ;I @C V;I ¼ rðDV;I rC V ;I Þ þ rT 2 @t K BT @C V ;I eq K VI ðC V C I C eq (1) V V C I Þ, @x K VI ¼ 4pðDV þ DI Þac expðDGIV =K B TÞ
(2)
where CV,I is actual concentration of vacancies and self-interstitials, DV,I is diffusivity of vacancies and self-interstitials, and KVI is the rate constant for recombination. The point defect distribution has been obtained using commercial software ‘‘FEMAG_CZ’’ [13]. To compare with actual experimental result, we used varying pull rate in our experiment. The concentration of point defects changed with varying pull rate, which affects micro defect formation and distribution. Initially, the thermodynamic property of point defect used the values of Sinno [14], and then these parameters were modified to match our experiment data. This step is needed because Tanahasi et al. [15] described that the interaction of point defects is greatly affected by impurity and thermal stress of the CZ process. The KLMC model has been widely used for thermal annealing and also used to describe the diffusion and clustering of vacancies in a Si lattice [16–18]. It is based on the Monte Carlo approach combined with the Poisson process [19]. The description of the diffusion and clustering of vacancies consists of many possible
Fig. 1. Multi-scale approach.
3593
events, which evolve as a series of independent events that are chosen according to their event rates in a Si system. The average rate of movement of defects in solids by thermal activation can be calculated by the classical rate theory. The rate of a migration event is then simply given by the Arrhenius expression as follows: E (3) ri ¼ n0 exp i , kB T where n0 is the attempting frequency as the prefactor and Ei is the migration energy, kB is the Boltzmann constant, and T is the simulation temperature. In this work we use 0.43 eV [9,10] for the migration energy of single vacancy jump. The bonding model used here to describe the vacancy cluster is the model developed by La Magna et al. [9]. Prasad and Sinno [20,21] mentioned that vacancies can interact up to the eighth nearest neighbor shell based on MD simulation results. The energy barrier equation for a vacancy hop is given by Dai et al. [22] as follows: 0 1 NN X jA @ DEi ¼ max 0; DEhop 0:5 DNBj Eb , (4) j¼1
where Ehop is the energy barrier for a single vacancy jump. NN is the nearest neighbor shell, and NBj is the change in the number of particle–particle bonds with interaction range j. The energy barrier is an important factor in KLMC simulation. The time interval is set by the inverse of the net hopping rate as follows: 1 , (5) q Rq P where Rq= i=0riq. In KLMC simulations, we can monitor the mean square displacement of the vacancies in the cell and evaluate the diffusivity dV as follows:
Dt ¼ P
dV ¼
N X 1 jRi ðtÞ Ri ð0Þj2 . lim t 6 t!1 i¼1
(6)
Here Ri(t) is the position of the ith atom or vacancy at time t. The point defect aggregation is well described by KLMC simulation. However, the initial condition of KLMC simulation is usually given by a general value. If we give more realistic condition for KLMC simulation, we could approach real experiment results. In the multi-scale method, we are using vacancy concentrations from the continuum model for initial vacancy concentrations of the KLMC model, so we could obtain the reasonable vacancy cluster distributions with the KLMC model.
3. Results and discussion We calculated point defect concentration with the continuum model. The simulations were performed under experimental conditions for 200 mm silicon crystal growth. The dynamic temperature was analyzed from 0 to 1000 mm for the length of crystal. CZ-Si crystal with the diameter of 200 mm was grown in the /1 0 0S direction with varying pull rates. The pulling rate is gradually decreased and then gradually increased between 0.9 and 0.6 mm/min. The obtained crystal was sliced off parallel to the vertical direction to observe the regions of grown-in defects. The defect regions were characterized by the lifetime mapping method and compared with the regions of the continuum model. To observe micro defect size and density, the obtained crystal was sliced off four wafers in the different pull rate region. Then, we analyzed the surfaces on wafers with MAGICS, which is an equipment to observe micro scale defects using the laser. The MAGICS detects defects above 20 nm well. The defect size and density were compared with the results of the atomistic model.
ARTICLE IN PRESS 3594
S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
Fig. 4. Vacancy diffusivity as a function of time for four temperatures.
Fig. 2. (a) Temperature profile and (b) point defect concentration (Cv–Ci) obtained from continuum model simulation with varying pull rates.
Fig. 3. Radial vacancy concentration estimated by the continuum model at 420, 570, and 720 mm along the crystal axis.
Fig. 2 shows the distributions for temperature and point defect concentration obtained from continuum model simulations with varying pull rates. Positive and negative values indicate the excess of self-interstitials and of vacancies, respectively. During the crystal growth process, the increase of pull rate leads to an increase of vacancy concentration. Due to the faster diffusion of self-interstitials than vacancies at high temperature, the selfinterstitial spreads out to the outer region of Si crystal. Fig. 3 shows the radial vacancy concentration estimated by the continuum model at three positions (420, 570, and 720 mm) along the length of the crystal obtained from continuum model simulation. These three positions coincide with the experimentally measured positions. The growth rate gradually increased from 0.60 to 0.90 mm/min. The vacancy concentration of 720 mm position is higher than that of 420 mm position
Fig. 5. (a) Number of clusters and (b) mean size of clusters as a function of time.
because of higher pull rate. The centers of 410, 570, and 720 mm positions are 7 1017, 1.5 1018, and 2.5 1018 cm–3, respectively. These values were used as the initial conditions of KLMC model simulations.
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
3595
In the KLMC model, long-range interactions were used to the eighth nearest neighbor, and the system temperature was clustering temperature range (1270–1370 K) [23,24]. The actual simulation volume is 2.0 10–17 cm3 and total simulation time is 0.01 s. The formation of vacancy clusters is well described by the KLMC model. Fig. 4 shows the diffusivity of vacancies as a function of time for variable temperatures such as the work by Haley et al.
[23]. The vacancy concentration in this case is 2.5 1018 cm3 at 720 mm position. The period of the initial constant diffusivity, indicating free vacancy diffusion, leads to decreasing diffusivity, as a power law of the form D(t)tr the exponent r=0.670.1, because clusters form. The diffusivities above 1470 K are almost constant, because clusters are not easily formed and most vacancies remain free. It means that a few vacancies agglomerate in clusters, but most are immediately dissociated to free vacancies. The diffusivity below 1370 K is decreased as the time increases. The diffusivity of vacancies remains greater at higher temperatures until 0.001 s, but it is decreased as the power law exponent is greater than the exponent at lower temperature after 0.001 s. Previous experiments [3,24,25] have shown that the size and the density of the vacancy defects are very sensitive to temperature. Especially, at the temperature of about 1370 K, large vacancy defects were increased and the numbers of small vacancy defects were decreased at the same time. At the temperature about 1270 K, small-size vacancy defects were more increased than large-size vacancy defects. Fig. 5 shows the formation of vacancy clusters for three different temperatures as a function of time. Fig. 5(a) and (b) shows the number of clusters and the mean size of clusters over time, respectively. At the temperature below
Fig. 6. Vacancy diffusivity as a function of time for three concentrations at (a) 1370, (b) 1320, and (c) 1270 K.
Fig. 7. (a) Number of clusters and (b) mean size of clusters as a function of time for vacancy concentrations.
ARTICLE IN PRESS 3596
S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
1370 K, the clusters that are formed at higher temperatures tend to be larger on average than those formed at lower temperature. The number of clusters increased and then gradually decreased as the time increased at the temperature about 1370 K, because small clusters are merged with other clusters or resolved into free vacancy. However, the number of clusters continuously increases as the time increases at the temperature of about 1270 K. Fig. 6(a)–(c) shows the diffusivity of vacancies for variable vacancy concentrations at the temperatures of 1370, 1320, and 1270 K, respectively. The Initial vacancy concentrations are the same as the vacancy concentrations obtained from the continuum model simulation at the center position of each crystal length. The diffusivity of vacancies decreases as clusters form, with a power law, where higher values correspond to higher vacancy concentration. Clusters formed at higher temperature are larger than clusters formed at lower temperature. So the diffusivity of vacancies decreases with a power law, where higher exponent values correspond to higher temperature. When the vacancy concentration is 0.7 1018 cm3, the diffusivity of vacancies at 1370 and 1320 K almost constantly remained. However, it was decreased as time increased at 1270 K. It means that the clustering needs the lower vacancy concentration in the lower temperature. Generally, smaller cluster has lower binding E, so the smaller clusters are easily resolved into free vacancy in the high temperature. It showed that the clustering phenomena is well affected by temperature and vacancy concentration. Fig. 7(a) and (b) shows the number of clusters and the mean size of clusters, respectively, at variable vacancy concentrations at 1270 K. The number of clusters increase continuously over time at 1270 K. It means that clusters are not easy to merge or resolve
Ingot position 420mm
Ingot position 570mm
Pixel Histogram
...40 ...20 ...10
41...
13
...40
61 66
...5 ...3
...10
54 67
...3
115 108
...1
Total = 643
...2 ...1
[PIzels]
Pixel Histogram
46 12 383
41... ...40 ...20
1586
37 63 2403 8440
...10
5485
...7 ...5
...2
[PIzels]
...20
119
...7
Ingot position 720mm
Pixel Histogram
40
41...
at low temperature. High vacancy concentration led to increase in the number of clusters. As shown in Fig. 7(b), while the mean size of clusters at 1270 K is slightly increased as time increases, it is irrespective of the vacancy concentrations. Fig. 8 shows defect distribution in the experiment estimated by MAGICS. Though we could not determine the exact size of micro defect from the data, pixel number roughly describes defect size. As the pull rate is increased, the total defect density increased with the density of the vacancy defects. Fig. 9(a) shows the comparison of the total defect density between experiments and simulation results. The higher the pull rate, the higher the vacancy concentration, so that it finally has a higher density of vacancy defects. Fig. 9(b) shows the radial cluster distribution obtained from the KLMC model simulations and experiments. To approach real experiment condition, we used the radial vacancy concentration of each crystal length obtained from the continuum model in Fig. 3. We could not directly present the actual defect size and density in KLMC results, because simulation volume is really different between simulation and experiment. For this reason, the total number of defects obtained from the KLMC simulations is different from that obtained from the experiments, but the corresponding defect distribution profiles are in good agreement with each other. Here, we only considered the interaction between the vacancies. Recent works for the defect formation of silicon have considered interaction between vacancies and impurities such as oxygen [26]. The interactions between the vacancies and impurities in silicon will be included in our further model, and their binding energies will be obtained from the ab initio calculations.
1540 1637
...5 ...3
4618 3852 Total = 19159
Fig. 8. Micro defect profile measured by MAGICS.
15847
...7
...2 ...1
[PIzels]
3166 3104 7062 5018 Total = 45140
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
3597
temperature region (1370–1270 K), larger clusters were generated at higher temperatures, and smaller clusters were generated at the lower temperatures. While the vacancy concentration led to increase in the number of clusters, the mean size of clusters was irrespective of vacancy concentrations. The larger clusters obtain higher binding energy, so clustering needs higher concentrations in higher temperatures, whereas it needs smaller concentrations in lower temperatures. Clustering phenomena were very susceptible to temperature and vacancy concentrations. The total number of vacancy clusters was linearly proportion to the crystal pull rate. The radial distribution of clusters obtained from multiscale simulations was in good agreement with the distribution of voids in experimental data. In this research, we tried to explain vacancy clustering phenomena, and we specifically researched vacancy behavior in the clustering temperature region (1370–1270 K). References
Fig. 9. Comparison between KLMC model simulations and experiments. (a) Total defect density. (b) Radial defect density.
4. Conclusion We have investigated vacancy cluster formation in silicon by means of the multi-scale method as an effective method for estimating micro defect in silicon crystal. The continuum model could effectively demonstrate the point defect distribution in real experiment conditions and describe the effects of the experimental parameters. KLMC simulations study for vacancy-assisted self-diffusion in silicon could provide important information on micro defect formation in silicon below melting temperature. At the temperature above 1470 K, the diffusivity was almost constant because the clusters were not easily formed. In the clustering
[1] R. Winkler, G. Behnke, in: H.R. Huff, W. Berghol, K. Sumino (Eds.), Semiconductor, ECS, Pennington, NJ, 1994, p. 973. [2] P.J. Rocksnoer, M.M.B. Van den Boon, J. Crystal Growth 53 (1981) 563. [3] V.V. Voronkov, J. Crystal Growth 59 (1982) 625. [4] R. Habu, T. Iwasaki, H. Harada, A. Tomiura, Jpn. J. Appl. Phys. 33 (1994) 1234. [5] R. Habu, A. Tomiura, Jpn. J. Appl. Phys. 35 (1996) 1. [6] T. Ebe, J. Crystal Growth 203 (1999) 387. [7] A. Bongiorno, L. Colombo, Phys. Rev. B 57 (1998) 8767. [8] F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. [9] A. La Magna, S. Coffa, L. Colombo, Nucl. Instrum. Methods B 148 (1999) 262. [10] L. Colombo, Physica B 273–274 (1999) 458. [11] R.A. Brown, et al., J. Crystal Growth 137 (1994) 12. [12] A. Voigt, C. Weichmann, J. Crystal Growth 266 (2004) 126. [13] N. Van den Bogaert, F. Dupret, J. Crystal Growth 171 (1997) 77. [14] T. Sinno, J. Crystal Growth 303 (2007) 5. [15] K. Tanahashi, M. Kikuchi, T. Higashino, N. Inoue, Y. Mizokawa, Physica B 273–274 (1999) 493. [16] L. Pelaz, L.A. Marques, M. Aboy, J. Barbolla, Appl. Phys. Lett. 82 (2003) 2038. [17] M. Aboy, L. Pelaz, L.A. Marques, L. Enriquez, J. Barbolla, J. Appl. Phys. 94 (2003) 1013. [18] R. Pinacho, P. Castrillo, M. Jaraiz, I. Martin-Bragado, J. Barbolla, J. Appl. Phys. 92 (2002) 1582. [19] K.A. Fichthorn, W.H. Weinberg, J. Chem. Phys. 95 (1991) 1090. [20] M. Prasad, T. Sinno, Phys. Rev. B. 68 (2003) 045206. [21] M. Prasad, T. Sinno, Phys. Rev. B. 68 (2003) 045207. [22] J. Dia, J.M. Kanter, S.S. Kapur, W.D. Seider, T. Sinno, Phys. Rev. B 72 (2007) 134102. [23] B.P. Haley, K.M. Beardmore, Phys. Rev. B 74 (2006) 045217. [24] B.M. Park, I.S. Choi, The ECS Proceedings Series, PV 2000-17, 2000, p. 19. [25] P.E. Blochl, E. Smargiassi, et al., Phys. Rev. Lett. 70 (1993) 2435. [26] J. Dia, W.D. Seider, T. Sinno, Mol. Simulat. 33 (2007) 733.
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Vacancy behavior in Czochralski silicon growth Sang Hun Lee a,b,, Jeong Won Kang c, Young Ho Hong a,d, Hyun Jung Oh a, Do Hyun Kim b a
Simulation Part, LG Siltron, Imsoodong, Gumi, Gyeongbuk 730-724, Republic of Korea Department of Chemical and Biomolecular Engineering, KAIST, Daejeon 305-701, Republic of Korea Department of Computer Engineering, Chungju National University, Chungju 380-702, Republic of Korea d Department of Ceramic Engineering Hanyang University, Seoul 133-791, Republic of Korea b c
a r t i c l e in f o
a b s t r a c t
Article history: Received 3 October 2008 Received in revised form 28 April 2009 Accepted 29 April 2009 Communicated by K.W. Benz Available online 14 May 2009
The vacancy defect behavior in silicon crystal growth has been investigated by kinetic Monte Carlo and continuum simulations. The vacancy concentration distributions in silicon crystal were obtained from continuum model simulations corresponding to experimental conditions and then the vacancy clusters distribution obtained from kinetic lattice Monte Carlo simulations. At the temperature above 1470 K, the diffusivity was almost constant because the clusters were not formed. In the clustering temperature region (1370–1270 K), the larger clusters were generated at the higher temperature, the smaller clusters were generated at the lower temperature. While the vacancy concentration led to increase in the number of clusters, the mean size of clusters was irrespective of the vacancy concentrations. Clustering phenomena were very susceptible to temperature and vacancy concentrations. The total number of vacancy clusters was linearly proportion to the crystal pull rate. The radial distribution of clusters obtained from multi-scale simulations was in good agreement with the distribution of voids in the experimental data. & 2009 Elsevier B.V. All rights reserved.
PACS: 61.50.Ah 61.72.Cc 61.72.jd Keywords: A1. Computer simulation A1. Diffusion A1. Kinetic lattice Monte Carlo A1. Point defects B2. Semiconducting silicon
1. Introduction The formation of defects in a silicon crystal during the growth process strongly affects the gate oxide integrity (GOI) characteristics of metal-oxide-semiconductor (MOS) devices [1]. Grown-in defects are formed by the aggregation of excess vacancies or selfinterstitials, and are thought to be related to the distribution of the point defect concentration during the crystal growth process. There are many theoretical and experimental studies for defect formation. Voids are formed by vacancy aggregation [2]. The type of grown-in defect depends on the growth rate. Voronkov [3] showed that the type of grown-in defect depends on the ratio V/G (where V is growth rate and G is axial temperature gradient of the crystal) and that there is a critical ratio that defines the defect type. Such studies prompted comprehensive modeling of point defect dynamics, and many researches have been developing and testing all scale simulation models including quantum mechanical calculations. Habu et al. [4,5] proposed models to describe the
Corresponding author at: Simulation Part, LG Siltron, Imsoodong, Gumi, Gyeongbuk 730-724, Republic of Korea. Tel.: +82 54 713 3799; fax: +82 54 713 3518. E-mail address: [email protected] (S.H. Lee).
0022-0248/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2009.04.044
dominant point defect species, and their model realistically describes the two-dimensional distribution of point defects in Czochralski (CZ) crystals. Ebe [6] performed analytical calculations on the behavior of point defects in the growing silicon crystal. Several theoretical descriptions of the vacancy cluster formation process in silicon exist. Bongiorno and Colombo [7] used the Stillinger–Weber potential [8] and molecular dynamics simulations to demonstrate that strain fields due to the distortion in the silicon lattice near vacancies affect the capture radius in different directions from a vacancy cluster. La Magna et al. [9] proposed a kinetic lattice Monte Carlo (KLMC) model for the description of vacancy diffusion and clustering phenomena. Colombo [10] described native defects and their interaction in silicon using several KLMC models. In this work, the interaction of the point defects in silicon crystal growth has been investigated by the multi-scale method combining KLMC and continuum models. Such a strategy is based on a coherent development of connected simulation tools with different accessibility domains. The continuum model can simulate whole processes but cannot consider atomic interactions directly whereas the KLMC model can simulate dynamic atomic interactions but the analyzing area is limited to nanometer scales. The multi-scale method is expected to play an important role in
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
the investigation of micro defects formation and distribution. In this work, we present the simulation results together with the experiment results corresponding to each simulation to validate our simulations.
2. Simulation methods Fig. 1 shows the process of multi-scale method, which consists of two types of simulation module for analyzing void generation with variable pull rate. Since the continuum model is useful for analyzing the whole system, point defect distributions of the whole crystal are obtained from continuum model simulation. However, this could not describe the aggregate of vacancies [11,12]. The KLMC model makes up for the weak point of the continuum model. This model is based on interactions between individual defects. These compliments are a good approach to understand void generation. Multi-scale method is a strategy to connect different scale simulation modules. We tried to analyze the real condition of industrial processes like changing pull rate. The continuum models for point defect dynamics have been developed by considering diffusion and recombination of defect species like vacancies, self-interstitials and impurities such as oxygen. We consider only native point defects in this continuum model. The concentrations of point defects during the Czochralski process are obtained by solving Eqs. (1) and (2) as follows: DV;I rC V;I Q V ;I @C V;I ¼ rðDV;I rC V ;I Þ þ rT 2 @t K BT @C V ;I eq K VI ðC V C I C eq (1) V V C I Þ, @x K VI ¼ 4pðDV þ DI Þac expðDGIV =K B TÞ
(2)
where CV,I is actual concentration of vacancies and self-interstitials, DV,I is diffusivity of vacancies and self-interstitials, and KVI is the rate constant for recombination. The point defect distribution has been obtained using commercial software ‘‘FEMAG_CZ’’ [13]. To compare with actual experimental result, we used varying pull rate in our experiment. The concentration of point defects changed with varying pull rate, which affects micro defect formation and distribution. Initially, the thermodynamic property of point defect used the values of Sinno [14], and then these parameters were modified to match our experiment data. This step is needed because Tanahasi et al. [15] described that the interaction of point defects is greatly affected by impurity and thermal stress of the CZ process. The KLMC model has been widely used for thermal annealing and also used to describe the diffusion and clustering of vacancies in a Si lattice [16–18]. It is based on the Monte Carlo approach combined with the Poisson process [19]. The description of the diffusion and clustering of vacancies consists of many possible
Fig. 1. Multi-scale approach.
3593
events, which evolve as a series of independent events that are chosen according to their event rates in a Si system. The average rate of movement of defects in solids by thermal activation can be calculated by the classical rate theory. The rate of a migration event is then simply given by the Arrhenius expression as follows: E (3) ri ¼ n0 exp i , kB T where n0 is the attempting frequency as the prefactor and Ei is the migration energy, kB is the Boltzmann constant, and T is the simulation temperature. In this work we use 0.43 eV [9,10] for the migration energy of single vacancy jump. The bonding model used here to describe the vacancy cluster is the model developed by La Magna et al. [9]. Prasad and Sinno [20,21] mentioned that vacancies can interact up to the eighth nearest neighbor shell based on MD simulation results. The energy barrier equation for a vacancy hop is given by Dai et al. [22] as follows: 0 1 NN X jA @ DEi ¼ max 0; DEhop 0:5 DNBj Eb , (4) j¼1
where Ehop is the energy barrier for a single vacancy jump. NN is the nearest neighbor shell, and NBj is the change in the number of particle–particle bonds with interaction range j. The energy barrier is an important factor in KLMC simulation. The time interval is set by the inverse of the net hopping rate as follows: 1 , (5) q Rq P where Rq= i=0riq. In KLMC simulations, we can monitor the mean square displacement of the vacancies in the cell and evaluate the diffusivity dV as follows:
Dt ¼ P
dV ¼
N X 1 jRi ðtÞ Ri ð0Þj2 . lim t 6 t!1 i¼1
(6)
Here Ri(t) is the position of the ith atom or vacancy at time t. The point defect aggregation is well described by KLMC simulation. However, the initial condition of KLMC simulation is usually given by a general value. If we give more realistic condition for KLMC simulation, we could approach real experiment results. In the multi-scale method, we are using vacancy concentrations from the continuum model for initial vacancy concentrations of the KLMC model, so we could obtain the reasonable vacancy cluster distributions with the KLMC model.
3. Results and discussion We calculated point defect concentration with the continuum model. The simulations were performed under experimental conditions for 200 mm silicon crystal growth. The dynamic temperature was analyzed from 0 to 1000 mm for the length of crystal. CZ-Si crystal with the diameter of 200 mm was grown in the /1 0 0S direction with varying pull rates. The pulling rate is gradually decreased and then gradually increased between 0.9 and 0.6 mm/min. The obtained crystal was sliced off parallel to the vertical direction to observe the regions of grown-in defects. The defect regions were characterized by the lifetime mapping method and compared with the regions of the continuum model. To observe micro defect size and density, the obtained crystal was sliced off four wafers in the different pull rate region. Then, we analyzed the surfaces on wafers with MAGICS, which is an equipment to observe micro scale defects using the laser. The MAGICS detects defects above 20 nm well. The defect size and density were compared with the results of the atomistic model.
ARTICLE IN PRESS 3594
S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
Fig. 4. Vacancy diffusivity as a function of time for four temperatures.
Fig. 2. (a) Temperature profile and (b) point defect concentration (Cv–Ci) obtained from continuum model simulation with varying pull rates.
Fig. 3. Radial vacancy concentration estimated by the continuum model at 420, 570, and 720 mm along the crystal axis.
Fig. 2 shows the distributions for temperature and point defect concentration obtained from continuum model simulations with varying pull rates. Positive and negative values indicate the excess of self-interstitials and of vacancies, respectively. During the crystal growth process, the increase of pull rate leads to an increase of vacancy concentration. Due to the faster diffusion of self-interstitials than vacancies at high temperature, the selfinterstitial spreads out to the outer region of Si crystal. Fig. 3 shows the radial vacancy concentration estimated by the continuum model at three positions (420, 570, and 720 mm) along the length of the crystal obtained from continuum model simulation. These three positions coincide with the experimentally measured positions. The growth rate gradually increased from 0.60 to 0.90 mm/min. The vacancy concentration of 720 mm position is higher than that of 420 mm position
Fig. 5. (a) Number of clusters and (b) mean size of clusters as a function of time.
because of higher pull rate. The centers of 410, 570, and 720 mm positions are 7 1017, 1.5 1018, and 2.5 1018 cm–3, respectively. These values were used as the initial conditions of KLMC model simulations.
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
3595
In the KLMC model, long-range interactions were used to the eighth nearest neighbor, and the system temperature was clustering temperature range (1270–1370 K) [23,24]. The actual simulation volume is 2.0 10–17 cm3 and total simulation time is 0.01 s. The formation of vacancy clusters is well described by the KLMC model. Fig. 4 shows the diffusivity of vacancies as a function of time for variable temperatures such as the work by Haley et al.
[23]. The vacancy concentration in this case is 2.5 1018 cm3 at 720 mm position. The period of the initial constant diffusivity, indicating free vacancy diffusion, leads to decreasing diffusivity, as a power law of the form D(t)tr the exponent r=0.670.1, because clusters form. The diffusivities above 1470 K are almost constant, because clusters are not easily formed and most vacancies remain free. It means that a few vacancies agglomerate in clusters, but most are immediately dissociated to free vacancies. The diffusivity below 1370 K is decreased as the time increases. The diffusivity of vacancies remains greater at higher temperatures until 0.001 s, but it is decreased as the power law exponent is greater than the exponent at lower temperature after 0.001 s. Previous experiments [3,24,25] have shown that the size and the density of the vacancy defects are very sensitive to temperature. Especially, at the temperature of about 1370 K, large vacancy defects were increased and the numbers of small vacancy defects were decreased at the same time. At the temperature about 1270 K, small-size vacancy defects were more increased than large-size vacancy defects. Fig. 5 shows the formation of vacancy clusters for three different temperatures as a function of time. Fig. 5(a) and (b) shows the number of clusters and the mean size of clusters over time, respectively. At the temperature below
Fig. 6. Vacancy diffusivity as a function of time for three concentrations at (a) 1370, (b) 1320, and (c) 1270 K.
Fig. 7. (a) Number of clusters and (b) mean size of clusters as a function of time for vacancy concentrations.
ARTICLE IN PRESS 3596
S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
1370 K, the clusters that are formed at higher temperatures tend to be larger on average than those formed at lower temperature. The number of clusters increased and then gradually decreased as the time increased at the temperature about 1370 K, because small clusters are merged with other clusters or resolved into free vacancy. However, the number of clusters continuously increases as the time increases at the temperature of about 1270 K. Fig. 6(a)–(c) shows the diffusivity of vacancies for variable vacancy concentrations at the temperatures of 1370, 1320, and 1270 K, respectively. The Initial vacancy concentrations are the same as the vacancy concentrations obtained from the continuum model simulation at the center position of each crystal length. The diffusivity of vacancies decreases as clusters form, with a power law, where higher values correspond to higher vacancy concentration. Clusters formed at higher temperature are larger than clusters formed at lower temperature. So the diffusivity of vacancies decreases with a power law, where higher exponent values correspond to higher temperature. When the vacancy concentration is 0.7 1018 cm3, the diffusivity of vacancies at 1370 and 1320 K almost constantly remained. However, it was decreased as time increased at 1270 K. It means that the clustering needs the lower vacancy concentration in the lower temperature. Generally, smaller cluster has lower binding E, so the smaller clusters are easily resolved into free vacancy in the high temperature. It showed that the clustering phenomena is well affected by temperature and vacancy concentration. Fig. 7(a) and (b) shows the number of clusters and the mean size of clusters, respectively, at variable vacancy concentrations at 1270 K. The number of clusters increase continuously over time at 1270 K. It means that clusters are not easy to merge or resolve
Ingot position 420mm
Ingot position 570mm
Pixel Histogram
...40 ...20 ...10
41...
13
...40
61 66
...5 ...3
...10
54 67
...3
115 108
...1
Total = 643
...2 ...1
[PIzels]
Pixel Histogram
46 12 383
41... ...40 ...20
1586
37 63 2403 8440
...10
5485
...7 ...5
...2
[PIzels]
...20
119
...7
Ingot position 720mm
Pixel Histogram
40
41...
at low temperature. High vacancy concentration led to increase in the number of clusters. As shown in Fig. 7(b), while the mean size of clusters at 1270 K is slightly increased as time increases, it is irrespective of the vacancy concentrations. Fig. 8 shows defect distribution in the experiment estimated by MAGICS. Though we could not determine the exact size of micro defect from the data, pixel number roughly describes defect size. As the pull rate is increased, the total defect density increased with the density of the vacancy defects. Fig. 9(a) shows the comparison of the total defect density between experiments and simulation results. The higher the pull rate, the higher the vacancy concentration, so that it finally has a higher density of vacancy defects. Fig. 9(b) shows the radial cluster distribution obtained from the KLMC model simulations and experiments. To approach real experiment condition, we used the radial vacancy concentration of each crystal length obtained from the continuum model in Fig. 3. We could not directly present the actual defect size and density in KLMC results, because simulation volume is really different between simulation and experiment. For this reason, the total number of defects obtained from the KLMC simulations is different from that obtained from the experiments, but the corresponding defect distribution profiles are in good agreement with each other. Here, we only considered the interaction between the vacancies. Recent works for the defect formation of silicon have considered interaction between vacancies and impurities such as oxygen [26]. The interactions between the vacancies and impurities in silicon will be included in our further model, and their binding energies will be obtained from the ab initio calculations.
1540 1637
...5 ...3
4618 3852 Total = 19159
Fig. 8. Micro defect profile measured by MAGICS.
15847
...7
...2 ...1
[PIzels]
3166 3104 7062 5018 Total = 45140
ARTICLE IN PRESS S.H. Lee et al. / Journal of Crystal Growth 311 (2009) 3592–3597
3597
temperature region (1370–1270 K), larger clusters were generated at higher temperatures, and smaller clusters were generated at the lower temperatures. While the vacancy concentration led to increase in the number of clusters, the mean size of clusters was irrespective of vacancy concentrations. The larger clusters obtain higher binding energy, so clustering needs higher concentrations in higher temperatures, whereas it needs smaller concentrations in lower temperatures. Clustering phenomena were very susceptible to temperature and vacancy concentrations. The total number of vacancy clusters was linearly proportion to the crystal pull rate. The radial distribution of clusters obtained from multiscale simulations was in good agreement with the distribution of voids in experimental data. In this research, we tried to explain vacancy clustering phenomena, and we specifically researched vacancy behavior in the clustering temperature region (1370–1270 K). References
Fig. 9. Comparison between KLMC model simulations and experiments. (a) Total defect density. (b) Radial defect density.
4. Conclusion We have investigated vacancy cluster formation in silicon by means of the multi-scale method as an effective method for estimating micro defect in silicon crystal. The continuum model could effectively demonstrate the point defect distribution in real experiment conditions and describe the effects of the experimental parameters. KLMC simulations study for vacancy-assisted self-diffusion in silicon could provide important information on micro defect formation in silicon below melting temperature. At the temperature above 1470 K, the diffusivity was almost constant because the clusters were not easily formed. In the clustering
[1] R. Winkler, G. Behnke, in: H.R. Huff, W. Berghol, K. Sumino (Eds.), Semiconductor, ECS, Pennington, NJ, 1994, p. 973. [2] P.J. Rocksnoer, M.M.B. Van den Boon, J. Crystal Growth 53 (1981) 563. [3] V.V. Voronkov, J. Crystal Growth 59 (1982) 625. [4] R. Habu, T. Iwasaki, H. Harada, A. Tomiura, Jpn. J. Appl. Phys. 33 (1994) 1234. [5] R. Habu, A. Tomiura, Jpn. J. Appl. Phys. 35 (1996) 1. [6] T. Ebe, J. Crystal Growth 203 (1999) 387. [7] A. Bongiorno, L. Colombo, Phys. Rev. B 57 (1998) 8767. [8] F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. [9] A. La Magna, S. Coffa, L. Colombo, Nucl. Instrum. Methods B 148 (1999) 262. [10] L. Colombo, Physica B 273–274 (1999) 458. [11] R.A. Brown, et al., J. Crystal Growth 137 (1994) 12. [12] A. Voigt, C. Weichmann, J. Crystal Growth 266 (2004) 126. [13] N. Van den Bogaert, F. Dupret, J. Crystal Growth 171 (1997) 77. [14] T. Sinno, J. Crystal Growth 303 (2007) 5. [15] K. Tanahashi, M. Kikuchi, T. Higashino, N. Inoue, Y. Mizokawa, Physica B 273–274 (1999) 493. [16] L. Pelaz, L.A. Marques, M. Aboy, J. Barbolla, Appl. Phys. Lett. 82 (2003) 2038. [17] M. Aboy, L. Pelaz, L.A. Marques, L. Enriquez, J. Barbolla, J. Appl. Phys. 94 (2003) 1013. [18] R. Pinacho, P. Castrillo, M. Jaraiz, I. Martin-Bragado, J. Barbolla, J. Appl. Phys. 92 (2002) 1582. [19] K.A. Fichthorn, W.H. Weinberg, J. Chem. Phys. 95 (1991) 1090. [20] M. Prasad, T. Sinno, Phys. Rev. B. 68 (2003) 045206. [21] M. Prasad, T. Sinno, Phys. Rev. B. 68 (2003) 045207. [22] J. Dia, J.M. Kanter, S.S. Kapur, W.D. Seider, T. Sinno, Phys. Rev. B 72 (2007) 134102. [23] B.P. Haley, K.M. Beardmore, Phys. Rev. B 74 (2006) 045217. [24] B.M. Park, I.S. Choi, The ECS Proceedings Series, PV 2000-17, 2000, p. 19. [25] P.E. Blochl, E. Smargiassi, et al., Phys. Rev. Lett. 70 (1993) 2435. [26] J. Dia, W.D. Seider, T. Sinno, Mol. Simulat. 33 (2007) 733.