Modeling of point defect formation in silicon monocrystals

Microelectronic Engineering 69 (2003) 641–645 www.elsevier.com / locate / mee

Modeling of point defect formation in silicon monocrystals V.A. Zabelin*, V.V. Kalaev Soft-Impact Ltd., 27 Engels Av., P.O. Box 33, St. Petersburg 194156, Russia

Abstract The prediction of characteristics of point defects in silicon crystals is the final goal of numerical simulation of the whole process of growth from the melt, recognized today as a helpful tool in technology optimization. We present an effective numerical model with an accurate description of defect incorporation and evolution in a growing crystal. The model is applied to a 300-mm diameter CZ silicon crystal with an investigation of the sensitivity of the predicted defect characteristics to input parameters.  2003 Elsevier B.V. All rights reserved. Keywords: Silicon; Point defect; Oxygen precipitate; Defect nucleation

1. Introduction High quality monocrystalline silicon is the basic semiconductor material for producing integrated circuits. In recent years, crystal homogeneity and microstructure properties have become of high importance, which demands the optimization of the whole technological process. Numerical simulation of crystal growth process [1–7] has become a fruitful tool for supporting the technology refinement. Using models of heat and mass transfer for a crystal growth system [1–4] one can study the influence of growth parameters on crystallization conditions and subsequent defect formation in silicon [4–8]. For silicon crystals produced by the Czochralski technique the primary imperfections are point defects: octahedral vacancy aggregates (so-called voids) and particles of amorphous oxygen precipitates SiO 2 , which are

*Corresponding author. E-mail address: [email protected] (V.A. Zabelin).

usually observed in CZ silicon due to a high oxygen concentration in the crystal. An optimal engineering model of defect incorporation and clusterization is still under discussion [4–11]. The present contribution is focused on the development of a numerical model coupling the computer efficiency with an accurate description of defect evolution in a growing crystal. The model contains novel aspects in the description of the defect nucleation rate and the oxygen precipitate kinetics.

2. Simulation of initial defect incorporation Initial lattice defects (single vacancies and selfinterstitial atoms) are the source for point defect clusterization during a crystal thermal treatment, and their concentrations determine the crystal quality to a large extent. The 2D model of initial defect incorporation into the crystal and their subsequent recombination in a hot region near the crystallization front could be presented in terms of the following equation with respective boundary conditions and parameters:

0167-9317 / 03 / $ – see front matter  2003 Elsevier B.V. All rights reserved. doi:10.1016 / S0167-9317(03)00357-5

642

V. A. Zabelin, V.V. Kalaev / Microelectronic Engineering 69 (2003) 641–645

≠C ]x 1V ? \Cx 5 \ ? (Dx \Cx ) ≠t 1 4p a r (Di 1 Dv )

S

D

DG eq 3 exp 2 ] sC eq i C v 2 Ci Cvd kT (1) where C eq is the equilibrium defect concentration x and Cx is the actual defect concentration, Dx is the defect diffusivity, V is the crystal pulling rate, a r is the recombination radius, and DG is the recombination free energy barrier (x 5 i, v for self-interstitials and vacancies, respectively). The values of thermophysical parameters of vacancies and self-interstitials are still under discussion (see, for example, Ref. [9]). In the present work, we used the constant set proposed in Ref. [8]. Strong variation of defect equilibrium concentrations and diffusivities by the reason of temperature distribution and strong nonlinear interdependence make this problem complicated and one should use high order of convective and spatial partial derivative operators approximation in the transport equations. A precise determination of thermal field and oxygen concentration in the silicon crystal is of high importance for the prediction of the defect crystal structure. 2D steady-state models of global heat transport and flows are conventional today [1–3]; however, only unsteady 3D modeling of turbulent melt convection seems to provide the valid data on

the geometry of the melt / crystal interface shape and oxygen concentration [10,12,13]. Fig. 1 gives results of the calculation of defect incorporation into the crystal and their subsequent recombination using the model described above. Two variants of input parameters were used in the calculations. First, the melt / crystal interface geometry was taken from 2D global heat transfer calculations using the same model as in Ref. [4], and, second, the interface geometry was obtained with the application of the 3D unsteady analysis [12,13]. Details about growth parameters of the 300-mm crystal under the consideration could be found in [13]. It is evident from the calculations (Fig. 1A and B) that the quality of input information is very important for the analysis of defect incorporation. The concentration of prevailing defects is really sensitive to the geometry of the melt / crystal interface. For the case with the interface taken from 3D analysis, the crystal is mostly vacancy rich without observing a significant interstitial rich zone as in the case of the interface obtained with using the input data of 2D steady modeling.

3. The model of point defect clusterization Certain characteristics of point defects formed during the crystal growth and subsequent wafer annealing are required for further making the integrated circuits (IC). Voids and oxygen precipitates near wafer surface can damage precise sub-micron

Fig. 1. Spatial distribution (A) and radial distribution (B) for A–A cross-section of initial defect concentration for 300 mm diameter crystals with interfaces obtained by the 2D and the 3D models. Pulling rate is 0.67 mm min 21 .

V. A. Zabelin, V.V. Kalaev / Microelectronic Engineering 69 (2003) 641–645

IC elements. Supporting the conventional experimental approach, there are a number of models of point defect nucleation and growth [4,7,8] based on the thermodynamic description of defect behavior. In this work we present a model of simultaneous formation and evolution of voids and oxygen precipitates, which allows the prediction of point defect concentrations and size distributions. The gain in the Gibbs free energy associated with point defect formation can be written as: GV 5 2 nfV 1 lV n 2 / 3 5 2 nkT lnsCv /C veqd 1 lV n 2 / 3 ] 5 2 4 / 3p R 3 rV fV 1 4p 2 / 3Œ3sV R 2

(2)

GP 5 2 nfP 1 lP n 2 / 3 2/3 5 2 nkT hg lnsCv /C veqd 1 lnsCo /C eq o dj 1 lP n ] 5 2 4 / 3p R 3 rP fP 1 4p 2 / 3Œ3sP R 2 (3)

for voids (2) and oxygen precipitates (3), respectively. These expressions are written in two notations: first, ‘n’ is the number of aggregated defects (vacancies for voids and oxygen atoms for precipitates, respectively); in the second notation, the gain is evaluated through the equivalent point defect radius. The first term in both expressions determines the contribution of volume free energy fk , associated with initial defect supersaturation (k 5V, P for voids and oxygen precipitates, respectively). The second term represents the contribution of the surface energy. The coefficients lk at the terms in the equations correspond to the actual interface energy for the k111l silicon plane sk in a way that volume of an equivalent sphere is equal to the volume of a octahedron containing ‘n’ elementary defects. Since an oxygen precipitate consists of amorphous SiO 2 , its oxygen atomic density rP differs from the silicon 22 23 atomic density rV 52310 cm . So, the simultaneous aggregation of vacancies and oxygen atoms is necessary to relax the crystal lattice strain. The molecular volume ratio h is determined as V(SiO 2 ) /V(Si) ¯ 2.2 and the proper coefficient g is equal (h 2 1) / 2. Taking into account void reactions as (Vn 1 v ↔ Vn 11 , Vn11 1 i →Vn ) and surface kinetics the void growth gV and dissolution dV rates can be expressed as:

S]Dd C D d (R) 5 4p R S] C d gV (R) 5 4p R V

2

v

2

v

643

D

Di V 2 ] C i,S , d Di V,eq V,eq C v,S 2 ] d i,S

V v,S

D

(4)

where C Vx,S is the surface concentration and C V,eq x,S is the surface equilibrium concentration of an initial defect. The reaction rates are determined by the surface defect mobility, and this dependence evaluates as the ratio Dx /d, where Dx is the defect diffusion coefficient and d is the jump distance, which is assumed to be of the order of the silicon lattice parameter. Combining surface kinetics with diffusive defect fluxes in the vicinity of void one can obtain in the steady spherical approximation the following expressions for both the growth and dissolution rates: 4p R 2 V,eq gV (R)5]]] hDvsR?C V,eq v,S 1d ?Cvd2DisR?C i,S 1d ?Cidj d (R1d ) 4p R 2 V,eq dV (R)5]]sDv C V,eq v,S 2Di C i,S d d

(5)

The total consumption of initial defects by voids can be evaluated as: 4p R 2 Dv V,eq ]]] Jv,V 5 C 2 C v,S d, R 1d s v 4p R 2 Di V,eq Ji,V 5 ]]]sCi 2 C i,S d R 1d

(6)

The description of oxygen precipitate evolution is similar to that of voids; but it is more complicated because of simultaneous aggregation of vacancies and oxygen atoms: Pn 1 o 1 g v↔Pn11 . Hence, a self-consistent problem of surface reaction kinetics and diffusive transport of both species should be considered. The precipitate growth rate gP and the dissolution rate d P could be written similar to expression (5) as: 4p R 2 Do P,eq gP (R) 5 ]]]sR ? C o,S 1 d ? Cod, d (R 1 d ) 4p R 2 Do P,eq d P (R) 5 ]]] C o,S d

(7)

The total initial defect fluxes to precipitates are evaluated similar to (6) as:

V. A. Zabelin, V.V. Kalaev / Microelectronic Engineering 69 (2003) 641–645

644

d k (R) ? (≠R / ≠n) and using an appropriate function R(n). After this substitution, one can derive a modified Fokker–Plank equation:

4p R 2 Dv P,eq Jv,P 5 ]]]sCv 2 C v,S d, R 1d 4p R 2 Do P,eq ]]] Jo,P 5 C 2 C o,S d R 1d s o

(8)

The equilibrium concentration of initial defects on the void or precipitate surface can be obtained by the extreme condition for the free energy of point defect ≠G k formation: ] ≠n 5 0. For oxygen precipitates, a connection between vacancy and oxygen concentrations can be derived as follows: P,eq

P,eq g

C o,S sC v,S

d

S

2lP eq eq 5 C o sC v dg exp ]]] 3n 1 / 3 kT

D

(9)

An additional equation could be obtained from the stoichiometry condition Jv,P 5 g Jo,P . So, the system of these non-linear equations determines the vacancy and oxygen surface equilibrium concentration. The formation of vacancy–oxygen complexes V O 2 playing a noticeable role in the stage of crystal cooling is also taken into account, following the approach presented in [11]. The evolution of voids and oxygen precipitates is governed by the Fokker–Plank equation for the size distribution function w k (n,t): ≠w k (n,t) ≠h k (n,t) ]]] 5 2 ]]], ≠t ≠n h k (n,t) 5 f gk (n(R)) 2 d k (n(R)) g w k 1 ≠ 2 ] ] hw k f gk (n(R)) 1 d k (n(R)) g j 2 ≠n

(10)

It is useful to convert this equation from the ‘n’space to the ‘R’-space by substituting Wk (R,t) 5 w k (n(R),t) ? (≠n / ≠R), j k (R) 5 gk (R) ? (≠R / ≠n), l k (R) 5

≠Wk (R,t) ≠Hk (R,t) ]]] 5 2 ]]], ≠t ≠R Hk (R,t) 5 f j k (R) 2 l k (R) g Wk 1 ≠ ] hWk f j k (R) 1 l k (R) g j 2 ]]] 2 8p R rk ≠R

(11)

The point defect nucleation rates, which are used as conditions for the Fokker–Planck equation, can be estimated by the Zeldovich theory of steady-state nucleation [14] and are expressed as: Ik 5 Dz (n * ) ZrV exps 2 G *k /kTd,

S S

2

1 d Gk Z 5 2 ]] ]] 2p kT dn 2 2 k

f 5 ]]] 12p kTG *k

D

1 ] 2

,

U D

1 ] 2

n 5n *

kT(dn / dt) Dz (n) 5 2 ]]] sdGk / dnd

(12)

where the values with asterisk correspond to critical nuclei parameters. Combining this expression with the equations for the loss rates of initial defects, which take into account recombination and consumption by point defect, one can obtain a complete system governing point defect formation and evolution. In Fig. 2, radial distributions of the concentrations of voids (A) and oxygen precipitates (B) are shown for the same crystal as in Fig. 1. A lower vacancy concentration for the case with the interface computed within the 2D model results to a higher point

Fig. 2. Radial distribution of point defect density for 300-mm diameter crystals with different interfaces.

V. A. Zabelin, V.V. Kalaev / Microelectronic Engineering 69 (2003) 641–645

defect concentration. The oxygen precipitate concentration especially increases near the border between the vacancy rich crystal core and the interstitial rich outer ring. In this zone, the vacancy supersaturation is insufficient for an active void formation, and the oxygen supersaturation can be the major driving force of point defect formation. This border corresponds to the well-known oxidation staking fault (OSF) ring of anomalous oxygen precipitation [7].

4. Conclusion We have presented an efficient numerical model for predicting defect characteristics in a growing silicon crystal, containing 2D model of initial defect incorporation and point defect evolution model which is based on several novel approaches in the description of the defect nucleation rate and the oxygen precipitate kinetics. The application of the model to 300-mm silicon CZ growth has shown that the results of the calculations are noticeably sensitive to input parameters, in particular, to the melt / crystal interface geometry, which indicates that models of heat and mass transport in growth systems should be of high accuracy.

645

References [1] E. Dornberger, E. Tomzig, A. Seidl, S. Schmitt, H.-J. Leister, Ch. Schmitt, G. Muller, J. Cryst. Growth 180 (1997) 461. [2] A. Lipchin, R.A. Brown, J. Cryst. Growth 216 (2000) 192. [3] V.V. Kalaev, I.Yu. Evstratov, Yu.N. Makarov, E.M. Smirnov, A.I. Zhmakin, in: ECCOMAS-2000, September 11–14, Barcelona, Spain, 2000, CD-ROM Proceeding, N677. [4] V.V. Kalaev, V.A. Zabelin, Yu.N. Makarov, Solid State Phenom. 82–84 (2002) 41. [5] T. Sinno, R.A. Brown, W. Ammon, E. Dornberger, J. Electrochem. Soc. 145 (1998) 302. [6] V.V. Voronkov, R. Falster, J. Appl. Phys. 86 (1999) 5975. [7] V.V. Voronkov, R. Falster, J. Cryst. Growth 194 (1998) 76. [8] T. Sinno, R.A. Brown, J. Electrochem. Soc. 146 (1999) 2300. [9] T. Sinno, Electrochem. Soc. Proc. 2 (2002) 212. [10] V.V. Kalaev, D.P. Lukanin, V.A. Zabelin, Yu.N. Makarov, J. Virbulis, E. Dornberger, W. von Ammon, Mater. Sci. Semicond. Process. 5 (2003) 369. [11] V.V. Voronkov, R. Falster, J. Cryst. Growth 204 (1999) 462. [12] I.Yu. Evstratov, V.V. Kalaev, V.N. Nabokov, A.I. Zhmakin, Yu.N. Makarov, A.G. Abramov, N.G. Ivanov, E.A. Rudinsky, E.M. Smirnov, S.A. Lowry, E. Dornberger, J. Virbulis, E. Tomzig, W. von Ammon, Microelectron. Eng. 56 (2001) 139. [13] V.V. Kalaev, D.P. Lukanin, V.A. Zabelin, Yu.N. Makarov, J. Virbulis, E. Dornberger, W. von Ammon, J. Cryst. Growth 250 (2003) 203. [14] Yu.I. Frenkel, Kineticheskaya Teoria Zhidkostey, Nauka, Leningrad, 1975, in Russian.