Numerical and experimental investigation of effect of oxygen concentration on grown-in defects in a Czochralski silicon single crystal

Grown-in defect-free wafers are required in silicon semiconductor devices. A point defect concentration simulation was performed along with an experimental investigation, demonstrating a wide range of oxygen concentrations from 1.6 × 1017 to 9.1 × 1017 cm−3 in crystals. Thus, the effect of oxygen atoms in a Czochralski silicon single crystal with grown-in defect behavior was revealed. Consequently, the increasing vacancy concentration trapped by the oxygen atom (oxygen coefficient) was estimated as 4.61 × 10−5 per oxygen atom. Previously, for obtaining the oxygen coefficient, a regression equation assuming thermal equilibrium concentrations of vacancy (V) and interstitial Si (I) was applied to the experimental results. However, the interface shape, thermal stress, and hot-zone structure of the experimental level needed to be arranged; this affected the grown-in defect behavior. In this study, the oxygen coefficient and thermal equilibrium concentration of V and I were determined uniquely without arranging the situations experimental level.


Introduction
In recent years, the demand for grown-in defect-free silicon (Si) wafers has increased for use in Si semiconductor devices with increasing integration and decreasing scale. Grown-in defects in Si crystals are formed by aggregating excess vacancies (V ) or self-interstitials (I) that are introduced at the solid-liquid interface. According to the Voronkov theory, 1,2) the defect type, i.e. V-or I-dominant, is determined by the v/G ratio, which is defined as the growth rate (v) over temperature gradient in the growth direction near the solidliquid interface (G). Thus, if the v/G value of a crystal is larger than the critical value, cri cri the defect type would be V-dominant. Conversely, if the v/G value is smaller than x , cri the defect type would be I-dominant. Hence, controlling the growth condition of v/G to x cri is necessary for obtaining a grown-in defect-free crystal.
Reportedly, x cri changes depending on the concentrations of impurity atoms such as dopants and oxygen incorporated in the crystal [3][4][5][6] and the shape of the solid-liquid interface during crystal growth. 7,8) Furthermore, it has been recently reported that x cri changes with thermal stress during crystal growth. [9][10][11][12][13][14][15][16][17][18] Therefore, to obtain the grown-in defect-free crystal, v/G must be adjusted to x cri during crystal growth according to each situation. Thus, it is important to understand the behavior of x cri with respect to impurity atoms, solid-liquid interface shape, and thermal stress.
In the impurity atoms, the change in x cri is caused by the excessive introduction of vacancies or interstitial Si trapped by the impurity atom in the crystal, and excess concentration is expressed as C trap ( a = C N, trap where C trap is the concentration of trapped V or I, N is the impurity concentration, and α is the impurity coefficient). 4) Therefore, various effects of impurity atoms on the α value must be realized.
There are two methods for determining the α value. In the first method, proposed by Voronkov et al., the α value can be obtained using a formula representing the relationship between x cri and N, that is, using x , cri,0 the critical v/G value without doping; x , cri the critical v/G value with doping; N; and the coefficient β; ,mp free 19,20) Here, C V,mp free and C I,mp free are the concentrations of free vacancy and free interstitial silicon at mp, respectively. The relationship between x cri and N can be obtained using the experimental results, while x cri,0 and β can be obtained via regression analysis. Finally, on the basis of the obtained β value and the assumed - ,mp free values, the α value can be determined. Using this method, the α value of each impurity can be obtained. 4) The second method entails the use of a point defect concentration simulation (PDSim), which can calculate grown-in defect distributions.
In the first method, to reveal the effect of impurity atoms on x cri and α, it is necessary to arrange solid-liquid interface shapes and thermal stress on the interfaces affecting x . cri However, it is difficult to control the oxygen concentration without changing the growth conditions and hot-zone (HZ) structures affecting the interface shapes. Therefore, the range of oxygen concentration used for investigation is extremely narrow.
The second method involves the use of PDSim. A wide range of oxygen concentrations can be investigated as changes in solid-liquid interface shapes and thermal stress can be considered to determine α. Although there are cases wherein PDSim was used to investigate the impact of α on the point defect behavior and tendency (x cri ), [21][22][23] there has been no investigation on determining α by using PDSim.
Therefore, in this study, to focus on oxygen atoms and investigate the effect of oxygen atoms on α, crystals with a wide range of oxygen concentrations were prepared, and α was determined using PDSim. Subsequently, α values determined by PDSim and those obtained using the regression equation were compared to examine the difference. Thus, we demonstrated the usefulness of using PDSim for calculating α.

Experiment
Czochralski silicon crystals with diameters of 200, 300, and 450 mm, different HZs, and oxygen concentrations from 1.6 × 10 17 to 9.1 × 10 17 cm −3 (conversion factor: 3.03 × 10 17 cm −3 , JEIDA), were prepared. The crystals were grown in the 〈100〉 direction, and the resistivity was 10-70 Ω•cm. Notably, it is well known that the change in resistivity (10-70 Ω•cm) has negligible effect on point defect behavior. 4) Each crystal was grown by gradually changing v to obtain the boundaries of grown-in defects. Consequently, 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. Post growth, the crystals were cut to a thickness of 1 mm along the growth direction. The P I −L/DL boundary (hereinafter referred to as the L/DL boundary) was obtained using the Cu decoration and Wright etch methods. 24) 2.2. Point defect concentration simulation method and determination of a and C V free 2.2.1. Point defect concentration simulation method. A point defect concentration simulation (PDSim), 24) which calculates the concentration distributions of V and I in the crystal during crystal growth, was used to investigate the effect of oxygen concentration on the grown-in defects behavior. Temperature distribution and thermal stress in the crystal were obtained considering the actual solid-liquid interface shape. 24) The mean stress P was used as the representative value of the stress owing to hydrostatic pressure and was defined as

rr zz
Here, a negative value of P indicates compressive stress, whereas a positive value indicates tensile stress.
The PDSim entails solving the advection-diffusion equation and considering the pair annihilation of V and I, as shown in Eq. (1). Next, the boundary condition is defined by assuming that the point defect concentrations at the solidliquid interface and crystal surface are thermal equilibrium concentrations: where C i is the actual point defect concentration of V or I, C i,eq is the thermal equilibrium concentration of V or I, J i is the diffusion flux of V or I, and K VI is the reaction coefficient of the pair annihilation of V and I. Moreover, under the existing temperature distribution in a silicon crystal, J and K VI are defined using Eqs. (2) and (3), respectively, as follows: where D i is the diffusion coefficient of V or I; Q i is the reduced heat transport of V or I, which is defined as the heat flux per unit flux of 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 a 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. 25) Furthermore, the thermal equilibrium concentration and the diffusion coefficient of V or I considering the thermal stress effect are given by Eqs. (4) and (5), respectively.
where C i 0, and D i 0, are the pre-factors, E i is the formation energy, and a i is the stress coefficient of V or I. The superscripts f and m indicate formation and migration, respectively.
The effect of oxygen is introduced into the PDSim, as shown in Eqs. (6) and (7): The total thermal equilibrium V concentration (C V total ) at the mp (T , mp solid-liquid interface) is the sum of free vacancy concentration (C V free ) and oxygen-trapped vacancy concentration (C V trap ), as shown in Eq. (6). Moreover, C V trap is assumed using the product oxygen coefficient (a) and oxygen concentration (C Oi ), as shown in Eq. (7). 4) Using the aforedescribed procedure for PDSim, the distributions of relative supersaturation at 1273 K (DC V ) were obtained for the crystal, as shown in Eq. (8). ,eq The value of DC V proposed by Tan et al. 26) and Brown et al. 27) was used to determine whether the point defects were V-or I-dominant and grown-in defect boundaries.
The physical properties used for PDSim are listed in Table I. Here, a f and a m were determined by the Suewaka and Nakamura, 24) whereas the other values were determined by Nakamura et al. 28) 2.2.2. Determination of a and C V free from PDSim. Experimental L/DL boundaries of crystal Nos. A-F were obtained, as shown in Sect. 2.1. As the L/DL boundary predicted by PDSim can be obtained by assuming the position of Table I. Physical properties used in PDSim. 3 24) the combinations of a and C V free consistent with experimental results were determined for each experiment.

Determination of a and C V
free from the regression equation Voronkov simplified the advection/diffusion and pair annihilation reactions of V and I and showed that x cri can be expressed as Eq. (9) using the physical properties of point defects. 1) Here, Eq. (9) is the case where the thermal stress effect is introduced.   Table II shows the specifications of the prepared crystals: crystal diameter, C , Oi interface height (h), and P. Here, h represents the difference between the position of the center and the outer periphery of the solid-liquid interface of the crystal. When h is positive, the crystal is concave, and when it is negative, the crystal is convex toward the melt side. Figure 1 shows a case of experimental results for crystal No. D. Here, (a) shows the relationship between the crystal position and the v/G value of the crystal center, and (b) shows the grown-in defect distribution in the crystal. In addition to the L/DL defect area, the Void, OSF, P V , and P I areas were obtained using the grown-in defects evaluation method. 24 decreased with increasing C , Oi as reported in previous studies, 3,4) and is accompanied by an increase in total vacancy concentration due to vacancies trapped by oxygen atoms, as shown in Eqs. (6) and (7).         Figures 5 and 6 show the results of the experiment and PDSim of crystals A and D, respectively. Here, (a) depicts the grown-in defect distribution obtained from the experiment, (b) shows the results of PDSim using the new physical properties, (c) shows the results of PDSim using the conventional physical properties, and (d) reflects the v/G value at the crystal center. In the PDSim results of (b) and (c), the iso-concentration line of DC V is drawn, and the blue line is the L/DL boundary predicted by PDSim.

Determination of α and
For crystal No. A, shown in Fig. 5, the L/DL boundary position predicted by PDSim using new physical properties shown in (b) agrees with the experimental results of (a), whereas the PDSim results using conventional physical properties shown in (c) disagreed with (a).
Crystal No. D's results (Fig. 6) showed that the L/DL boundary positions of both (b) and (c) were in agreement with (a).
Thus, L/DL boundary positions predicted by PDSim, using the new physical properties, were in agreement with the experimental results for all crystals. However, PDSim results used conventional physical properties depending on the crystals (oxygen concentrations). In summary, Fig. 7 shows the relationship between each oxygen concentration (crystal) and the ratio of the difference between   When new physical properties were used, it showed good agreement with the experimental results within ±1% regardless of the oxygen concentration. However, when the conventional physical properties were used, the error was approximately 5% for high oxygen concentrations ( >´-C 8 10 cm 3 ). At < <´-C 5 10 7 10 cm , 17 Oi 17 3 the error was nearly 0%. This is because the conventional physical properties of ( ) were obtained by fixing α = 6.35 × 10 −5 , and the investigation range of C Oi was as narrow as 6 × 10 17 -8 × 10 17 cm −3 . 24) Therefore, the results of PDSim with new physical properties, obtained in this study, showed a good agreement with the experimental results over a wide range of oxygen concentrations. In particular, the agreement of experimental and PDSim results with low oxygen concentrations can have a significant impact on grown-in defect control for Si power devices. 29,30) 3.2.1.2. Application in other experimental results. New physical properties were applied to the experimental results of Nakamura et al. 4) and verified for universality. However, our estimation method of G used to calculate x / L DL exp.
is different using the method in this study. Therefore, G in the experimental results was recalculated using our method.  However, using the recalculated result of G, a =´-4.31 10 5 was obtained. This value was in agreement with that obtained using PDSim in this study, as shown in Table III. Red circles denote the results of new physical properties, and black squares denote the results of conventional physical properties.   indicates the values when P is minimum(crystal No. A) and maximum(crystal No. F), respectively. As shown in Table V, α values could obtained as 5.81 × 10 −5 and 6.31 × 10 −5 , respectively, and α values changed with stress. While using Eq. (10), D ( ) C T P , mp between the experimental levels must be constant; thus, the P value must be constant across all experiments. Therefore, the experimental results from Nakamura et al. 4) can reasonably confirm the oxygen concentration effect on x / L DL as the interface height and HZ are constant regardless of P.
Evidently, our estimation method for α using PDSim is useful because it can be obtained without fixing the growth conditions such as the interface height, P, and HZ structure.

Conclusions
The point defect concentration simulation for silicon crystal growth was performed to reveal the impact of an oxygen atom on the point defect behavior. Consequently, an oxygen coefficient, α = 4.61 × 10 −5 (increase factor in vacancy concentration per oxygen concentration) was determined uniquely without arranging the interface shape, thermal stress, and HZ structure at the experimental level. Additionally, a wide range of oxygen concentrations was investigated. The method using PDSim is useful as it enables the investigation of the effect of oxygen as well as other impurity atoms on the point defect behavior in wide concentration ranges without adjusting the experimental conditions.