Effect of thermal stress on point defect behavior during single crystal Si growth

Silicon devices currently require silicon wafers that are free of grown-in defects. Point-defect simulation for silicon crystal growth, which is performed using the advection-diffusion equation considering the pair annihilation of vacancies and self-interstitials, is improved in this study by considering the effect of thermal stress during the growth process. This effect is introduced into the point-defect simulation as the stress term of the formation enthalpy. The stress coefficients in the stress term are determined by analyzing the correlation between the grown-in defect patterns obtained from the actual pulling tests and simulations. The defect patterns obtained from the simulations performed using the improved point-defect simulation method in this study were congruent with the experimental results, particularly for crystals with large diameters and high stress. The improved simulation method proposed in this study shall be useful for designing new hot-zone conditions and enhancing the understanding of point-defect behavior.


Introduction
At present, grown-in defect free silicon wafers are required for silicon devices with decreasing scale and high integration. Grown-in defects in silicon crystals are formed by the aggregation of excess vacancies (V ) or self-interstitials (I).
According to the Voronkov theory, 1,2) whether the defect type is V-or I-dominant is determined by the v/G ratio, defined as the growth rate (v) over the growth direction temperature gradient near the crystal-melt interface (G). In other words, a crystal that has a v/G larger than the critical value, (v/G) cri = ξ cri , will become V-dominant, whereas those with v/G smaller than ξ cri will become I-dominant. Moreover, to obtain a grown-in defect free crystal, it is necessary to control the growth condition of v/G to ξ cri . 3) On the basis of Voronkov theory, Vanhellemont claimed the significance of the effect of thermal stress near the crystal-melt interface on ξ cri . [4][5][6][7] Recently, the impact of thermal stress on ξ cri was determined in detail by Sueoka et al. using density functional theory (DFT): the compressive stress in the crystal near the crystal-melt interface increases the V concentration, thus resulting in a decrease in ξ cri . [8][9][10][11][12] Nakamura et al. also reported experimental evidence that compressive stress resulted in a shift to V-dominance and a decrease in ξ cri . 13) Sueoka et al. and Kamiyama et al. introduced the thermal stress effect on the point defect simulation, wherein the concentration distributions of V and I in the crystal during the growth process were calculated by solving the advection diffusion equation and considering the pair annihilation of V and I. 14,15) On the other hand, Abe et al. put a unique interpretation on the determination of defect type V-or I-dominant from the experimental data, i.e. the growth interface is filled with only vacancies, and self-interstitials are generated when the thermal stress exceeds a critical value. In other words, they claimed that the compressive stress results in a shift to I-dominant. 16,17) However, in these reports of Sueoka et al. [8][9][10][11][12] and Nakamura et al. 13) the point defect behavior due to thermal stress was confirmed by the change in ξ cri . In the studies by both Kamiyama et al. 14) and Sueoka et al. 15) it was difficult to verify the impact of stress effect on point defect behavior. This was because of the low stress conditions and no comparison between experimental and simulation grown-in defect patterns in Sueoka's study 15) and due to the discrepancy between the actual crystal-melt interface shape and the simulation in Kamiyama's study. 14) Similarly, in the report of Abe et al. 16,17) verification of the assumption is difficult as the actual and critical stresses were not shown. Therefore, in this study, the effect of thermal stress is introduced on the point defect simulation, and its influence on the grown-in defect behavior is investigated by comparing experimental results and simulation ones. In addition, the physical properties required to introduce the thermal stress effect on a point defect simulation have been determined.

Physical model
There are two parts in this section: the first subsection describes a well-known point defect simulation method, whereas the second introduces the effect of thermal stress on the point defect simulation.

Point defect simulation
To obtain accurate results for the grown-in defect pattern, the point simulation method used in this study involves three steps described below.
2.1.1. Prediction of the temperature distribution in the crystal. The temperature distribution of the crystal, particularly near the crystal-melt interface, has a significant impact on the point defect behavior, thus an accurate estimation of the temperature distribution is extremely important. Generally, temperature distributions can be calculated using a global heat transfer simulator, which are commercial software programs such as CGSim, 18,19) FEMAG, 20) and CrysMAS. 21) However, the shapes of crystal-melt interface obtained by these simulators differ slightly from actual situations; thus, the accuracy of temperature distributions obtained with this method is insufficient for use in point defect simulation. Therefore, to estimate the temperature distribution near the crystal-melt interface more accurately, the temperature distributions are calculated using the following method. First, the reference temperature distribution on the crystal surface is obtained using the global heat transfer simulator CGSim. Second, the actual crystal-melt interface shape is measured based on the growth striations in an actual crystal. Finally, the reference temperature on the crystal surface and the crystal-melt interface shape with a melting point of 1685 K are defined as the boundary conditions. Then, temperature distributions are recalculated by considering the thermal advection accompanying crystal growth. 2.1.2. Prediction of thermal stress. The thermal stresses (σ rr , σ θθ , σ zz ) in the crystal are estimated with ABAQUS using the temperature distribution obtained through the method described in Sect. 2.1.1. Three elastic constants (C11, C12, and C44) with temperature dependence 22) are converted to Young's modulus, shear modulus, and Poisson's ratio, respectively, for an elastic isotropic body using the Reuss average method. 22) The coefficient of thermal expansion reported by Okada et al. 23) is used. The mean stress P is used as a representative value of the stress due to hydrostatic pressure and is defined as rr zz ( ) / Here, a negative value of P indicates a compressive stress, whereas a positive value indicates a tensile stress.  1) and (2). [24][25][26] Then, the boundary condition is defined by assuming that the point defect concentrations at the crystal-melt interface and crystal surface are thermal equilibrium concentrations: where C V and C I are the concentrations of V and I, respectively; J V and J I are the diffusion fluxes of V and I, respectively; v is the growth rate (moving speed); K VI is the reaction coefficient of the pair annihilation; and C V eq and C I eq are the thermal equilibrium concentrations of V and I, respectively. The first, second, and third terms on the right side of these equations describe diffusion, advection, and pair annihilation rate, respectively. The parameters J and K VI are explained in greater detail below. Generally, under the existing temperature distribution in a silicon crystal, J is defined for V and I by Eqs. (3) and (4), respectively, whereas K VI is defined by Eq. (5) where D V and D I are the diffusion coefficients of V and I, respectively; Q V and Q I are the reduced heat transport of V and I, respectively, which is defined as the heat flux per unit flux of a component atom in the absence of a temperature gradient; T is the absolute temperature; k B is the Boltzmann constant; a c is the critical distance (=0.543 nm) within which pair annihilation reaction can occur between V and I; and DG IV is the barrier energy for pair annihilation, as obtained by Sinno et al. 27) The second term on the right side of Eqs. (3) and (4) represents the effect of temperature gradient on diffusion.

Introducing the thermal stress effect on the point defect simulation
The thermal equilibrium concentrations and diffusion coefficients are given by Eqs. (6)- (9). From these, it can be implied that the thermal stress in the crystal changes both the thermal equilibrium concentration and the diffusion coefficient, because the formation enthalpy, H , f and migration enthalpy, H , m both contain a stress term (aP). In the conventional point defect simulation method, H uses the formation energy, E, whereas ignoring aP, as the value of aP is much smaller than E. However, it has been reported that aP has a significant impact on point defect behavior, 4,[7][8][9][10][11][12] and thus the simulation method presented here considers the thermal stress effect on the enthalpies of formation and migration = - where C 0 and D 0 are the pre-factors; a is the stress coefficient; and the superscripts f and m indicate formation and migration, respectively.  Table I. 28) As new unknowns (a ,  9) to describe the thermal stress effect, these values must be first determined.
The Voronkov criterion x P cri ( ) 1,2) used to introduce the thermal stress effect can be written as follows: where T mp is the melting point (1685 K). A change in x P cri ( ) indicates a change in the point defect behavior induced by the stress effect, and thus the stress effect on the point defect behavior can be predicted by comparing the difference between x P cri ( ) and x 0 . cri ( ) Here, the value of x 0 cri ( ) can be calculated as 0.163 mm 2 min −1 K −1 based on the physical properties given in Table I.
However, there are no unique solutions for the four unknown properties (a , ) the stress effect on the migration enthalpy can be ignored, i.e. it can be assumed that 9) it can be calculated that x -10 cri ( ) = 0.153 mm 2 min −1 K −1 . Therefore, x -10 cri ( ) decreases by 6.1% compared with x = 0 0.163 . cri ( )( ) It is necessary to control the actual value of v/G within ±4% for the manufacture of grown-in defect free crystals, and thus the change in the critical v/G due to the stress effect is greater than the actual v/G change to be controlled. Furthermore, the influence of the combination of a V f and a I f on the change in x -10 cri ( ) is investigated. Figure 1 shows the relationship between x -10 cri ( ) and and 5, respectively) were grown with gradually decreasing growth rates to obtain grown-in defect patterns from V-to Idominant. All of the experiments were performed with different HZ types. As a result, grown-in defects were obtained from voids, oxidation-induced stacking faults (OSF), P V , P I , and L/DL (dislocation clusters). Here, P V and P I are the defect-free regions with P V being V-dominant and P I being I-dominant. For the crystals grown, the oxygen concentration was 10-12 × 10 17 cm −3 (Old ASTM, F121-79, 1979) and the resistivity was 10-30 Ω cm.
Samples were cut along the growth direction to evaluate the grown-in defect patterns. The defect boundaries were obtained using the minority carrier life time-map and visual and microscopic observations after the application of heat treatment and/or Cu-decoration and selective etching, as summarized in Table II. 3.1.3. Parameter fitting method. The fitting method is shown in Fig. 2; in this study, the P I -L/DL boundary is employed. Figure 2 Table I  Comparing the experimental result in Fig. 2(a) with the simulation result in Fig. 2(b) reveals that there is no a V f value that can concurrently satisfy the P I -L/DL boundary of the crystal center position and its radial profile because the C V 0, given in Table I was obtained without introducing the stress effect.
Thus, C V 0, must be adjusted for each a V f so that the P I -L/DL boundary of the crystal center matches between the  experimental and the simulation results, as shown in Fig. 2(c). In addition, for each combination of a V f and C , V 0, the best matching condition that agreed with the experimental results for both the P I -L/DL boundary of the crystal center position and its radial profile was determined.
In addition, the effect of oxygen atoms on a thermal equilibrium concentration of V was reported. 31,32) C V 0, was determined by considering the effect of oxygen atoms; hence, the thermal equilibrium concentration of V is obtained by rewriting Eq. (6) as follows: where C V i 0, ,O is the pre-factor introducing the effect of oxygen atoms.
The change in the thermal equilibrium concentration of V due to the oxygen atoms at T mp is given by Eq. (12): 31,32) = )is the thermal equilibrium concentration of V when the oxygen concentration is 0; Oi is the oxygen concentration defined by the old ASTM F121-79; and a is the increase in vacancy concentration relative to the oxygen concentration (a = 4 × 10 −12 ). 31,32) From Eqs. (6′) and (12), (13) and (13′) are derived. As a result, C V 0, with consideration on the effect of oxygen atoms was obtained These operations were performed for five experiments, and the combination a V f and C V 0, was obtained. As seen in Fig. 3(e), the P I -L/DL boundary patterns obtained from the simulation were not affected by changes in Table II. Evaluation flow for the grown-in defect patterns.  Table I; (c) simulation result introducing the stress effect by changing a V f using the physical properties in Table I and adjusting C .

Detection method OSF region
the effect of the thermal stress is not observed. This can be attributed to the uniform radial profile of the thermal stress as shown in Fig. 3(c). Its low value is due to the small crystal diameter (200 mm).
On the other hand, the cases for the 300 and 450 mm diameter crystals are different, as shown in Figs. 4 and 5. From Figs. 4(c) and 5(c), the radial thermal stress profiles were non-uniform, and their values are high at the center. These features are due to large crystal diameters and the concave interface shape in Figs. 4(a) and 5(a). Thus, the P I -L/DL boundary patterns, obtained from the simulation, change with changes in a , V f as shown in Figs. 4(e) and 5(e). Therefore, the large diameter of these crystals (300 and 450 mm) has a marked impact on the point defect behavior. In the case of Fig. 4, a V f = 0.11 meV MPa −1 and = C V 0, -3.65 10 cm 26 3 were obtained as the best fitting parameters. Table III lists 13) changes in ξ cri by 1 MPa of compressive stress obtained the relationship between ξ cri and thermal stress showed around 1.1 × 10 −3 mm 2 min −1 K −1 . Changes of ξ cri were estimated to be 1.0 × 10 −3 and 1.6 × 10 −3 from DFT analysis, when assuming the isotropic stress and plane stress respectively. 10) By employing the point defect simulation method in this study, we found out that stress coefficients a a and

Verification of the stress effect
In this subsection, the influence of thermal stress on the grown-in defect pattern is discussed. Figures 6 and 7 show the results of tests No. 1 and 3, respectively. In Figs. 6 and 7, (a) shows the experimental result of the grown-in defect pattern; (b) shows the simulation result obtained using the conventional point defect simulation method by applying the physical properties in Table I with no consideration of the thermal stress effect; and (c) shows the simulation result introducing the thermal stress effect and using the physical properties in Table I, Table I with adjusted C .
V,0 shifted, and the changes in the defect boundary shape were small. This is because the radial stress profile was uniform and its value was low, as shown in Fig. 3(c); and C V 0, was not adjusted in Fig. 6(b). As a result, the simulation results with and without introduction of the stress effect were almost the same, and these results are consistent with the experimental result shown in Fig. 6(a). Therefore, the defect pattern obtained using the conventional simulation method without introduction of the stress effect was in good agreement with the actual defect pattern for a crystal with a diameter of <200 mm and low stress. This result likely demonstrates the reason that the conventional point defect simulation method has produced results that are in good agreement with the experimental patterns. 28) Figure 7 shows the results for test No. 3 with a 300 mm diameter crystal and high stress; Fig. 7(d) shows the radial  Table I with adjusted C .  Table I with  compressive stress in the radial profile, and the compressive stress causes an increase in the thermal equilibrium concentration of V due to the decrease in the formation enthalpy as described in Eq. (6′). Thus, in this case, ΔC V on the crystal center side in Fig. 7(c) shifts to V-dominant compared to Fig. 7(b). As a result, the grown-in defect pattern obtained using the point defect simulation and considering the stress effect was in good agreement with the experimental result. Therefore, these results shows that the compressive stress leads to V-dominant. It should be noted that the introduction of the stress effect into the point defect simulation is extremely significant, particularly when there is high thermal stress in the crystals.

Conclusions
The point defect simulation for silicon crystal growth, which is calculated using the advection diffusion equation considering the pair annihilation of V and I, was improved in this study by considering the effect of thermal stress during the growth process. The effect of the thermal stress was introduced into the point defect simulation as the stress term of the formation enthalpy. The stress coefficients were deduced from the relationship obtained through the correlation of the grown-in defect patterns in actual pulling tests and the simulations:  Using the improved point defect simulation method in this study, the defect patterns obtained from the simulation were in good agreement with the experimental results. Thus, the improved simulation method presented here will be useful for designing new hot zone conditions and increasing the understanding of point defect behavior.  Table I and (c) the improved simulation method introducing the stress effect and using the physical properties in Table I Table I; (c) the improved simulation method introducing the stress effect and using the physical properties in Table I