- Email: [email protected]
Engineering analysis of microdefect formation during silicon crystal growth
Journal of Crystal Growth 225 (2001) 97–109
Engineering analysis of microdefect formation during silicon crystal growth Robert A. Brown*, Zhihong Wang, Tatsuo Mori Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave 3-208, Cambridge, MA 02139-4307, USA
Abstract The quality of single crystalline silicon wafers used as substrates for microelectronic devices is measured in terms of the type, size and number density of microdefects formed during crystal growth and subsequent processing. Native point defects}vacancies and self-interstitials}and impurities, such as oxygen and carbon, in the silicon crystal convect, diffuse, react and aggregate, driven by species super-saturation, to form defect structures at densities that are microscopically visible. We model defect dynamics using diffusion–reaction theory in which defect types are represented as chemical species and have concentrations that are governed by conservation laws, chemical thermodynamics and kinetics, using self-consistent values of equilibrium, transport and kinetic parameters for point defects and impurities. We refer to this approach as crystal defect dynamics analysis. The purpose of this paper is to demonstrate that this approach leads to a wide range of predictions that are consistent with experiments.The results presented here focus on the two-dimensional predictions of microdefect structures under growth conditions where both voids (cluster of vacancies) and self-interstitial aggregates are seen in regions of the crystal separated by the oxidation-induced stacking fault ring (OSF-ring). The location and structure of the OSF-ring is predicted in terms of point defect dynamics near the melt/crystal interface, which results in an annular region with almost balanced point defect concentrations (centered at r RD¼0 ), and in terms of a peak (centered at r RCV;max ) in the residual vacancy concentration after cluster formation. This residual vacancy concentration facilitates oxide precipitation during crystal growth, which seeds stacking fault formation during subsequent oxidation. Calculations show that RCV;max 4RD¼0 with both values scaling with V=GðRÞ, where V is the crystal pull rate and GðRÞ is the axial temperature gradient at r ¼ R along the melt/crystal interface. An annulus just outside the OSF-ring (centered at r Rfree ; RCV;max 4Rfree 4RD¼0 ) is identified where almost no microdefects are present, thereby identifying the temperature field that leads to almost perfect silicon. # 2001 Elsevier Science B.V. All rights reserved. Keywords: A1. Computer simulation; A1. Defects; B2. Semiconducting silicon
1. Introduction In the last decade there has been tremendous progress in the understanding of the cause for the *Corresponding author. Tel.: +1-617-253-5726; fax: +1617-253-8549. E-mail address: [email protected] (R.A. Brown).
myriad of crystallographic imperfections that are present in the otherwise perfect single-crystal silicon substrates grown by the Czochralski (CZ) crystal growth method. Today it is becoming well understood that the microdefects observed in asgrown silicon boules and the precipitation of oxygen and other impurities during growth and wafer processing is controlled by the dynamics of
0022-0248/01/$ - see front matter # 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 1 ) 0 0 8 2 5 - 9
98
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
intrinsic point defects}vacancies and self-interstitials}which diffuse, react and aggregate to form the structures observed by etching methods [1,2], particle counting [3], laser scattering [4] and other techniques. A new approach has been developed with the aim of connecting the dynamics of point defects and microdefect evolution with processing conditions during crystal growth and subsequent wafer processing; we call this approach crystal defect dynamics analysis. The role of intrinsic point defects in microdefect formation has long been postulated in the formation of microdefects such as the self-interstitial aggregates known as A- and B-defects [5,6] and the vacancy dominated D-defects in float zone silicon [7]. The transition between vacancy-rich defects at the core of a Czochralski crystal and selfinterstitial dominated defects near the periphery is shown in the X-ray topograph of the axial section of a CZ crystal in Fig. 1, taken from [8]. The experiment shows the change in the defect structure with variation in pull rate. The sharp band separating voids and stacking faults is denoted and responds to pull rate changes, moving inward for decreasing crystal pull rate and outward with increasing pull rate.
Fig. 1. X-ray topograph of crystallographic defects in axial cross-section of CZ crystal. The crystal was grown with variations in the pull rate, as shown by the scale on the left, The migration of the OSF-ring with pull rate is clearly seen.
Voronkov [9] was the first to connect theoretically point defect dynamics and crystal growth operating conditions, such as crystal pull rate V and temperature field, with observable transitions in microdefect structures during crystal growth. His approach has pioneered comprehensive modeling of point defect dynamics [8,10–13], which is the subject of this paper. The analysis used in this approach is built on the hypothesis that microdefects in silicon crystals are the result of the transport, equilibrium and reaction of intrinsic point defects and impurities. The reactions may be to either annihilate defects, as in the case of selfinterstitial and vacancy recombination, or to nucleate and aggregate defects to form larger, more complex species. Large clusters form into microscopically observable defects driven by large point defect (and impurity) super-saturations that result as the crystal cools and during wafer processing. Subtle changes in the temperature field during processing [8] and changes in point defect dynamics caused by interactions with impurities, such as boron and carbon [14–16] swing the balance from vacancies to self-interstitials and change the final microdefect structure. The success of point defect dynamics for the prediction of defect structures in silicon has been phenomenal. The predictions include understanding the sensitivity to the crystal temperature field and pull rate of important features of the microdefects in the crystal. These features include; the boundary between the vacancy- and selfinterstitial-rich regions in the crystal, the nucleation temperature, the defect number density and average observable size of microvoids in vacancydominated crystal growth and the onset temperature and density of self-interstitial defects. Most recently, the predictions have been extended to many of the features of oxygen precipitation observed during crystal growth and wafer annealing [17]. Most remarkably, a self-consistent set of predictions that are in reasonable agreement with experiment can be obtained using thermophysical properties for point defect equilibrium and transport grounded in atomistic simulations and fit to experimental data for the dynamics of the oxidation-induced stacking fault ring (OSF-ring), as described elsewhere [8,13].
99
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
The defect dynamics models used to describe these processes are built on simple concepts for transport and diffusion-controlled reaction for dilute species in a crystalline solid. The formulation is made more complicated by the need to describe clusters of point defects and impurities that range in size from dimers, trimers, etc. to clusters as large as 1400 nm in diameter, which contain Oð1010 Þ point defects! This range is achieved by coupling the discrete rate equations that describe the reactions of small clusters with the Fokker–Planck or continuous rate description for the dynamics of large species [11]. The formulation of the defect dynamics models is described in detail elsewhere and is only reviewed in Section 2 [11–13]. These formulations lead to extremely complex combinations of diffusion– reaction equations for the defect species representing monomeric species and aggregates. These equations are similar to equations describing classical chemical reactions [18]. We have developed novel numerical methods for the accurate and efficient parallel solution of these equations for both steady and transient dynamics [19]. Solutions and experimental evidence have suggested asymptotic analysis for specific situations that has provided valuable insight into scaling [8]. The purpose of this paper is to demonstrate the application of defect dynamics analysis for several specific problems that arise in Czochralski crystal growth. Of specific interest are the details of the formation of the OSF-ring in large-scale crystals. We use the analysis of this microdefect structure as the basis for the presentation that follows. The experimental features of the OSF-ring are reviewed in Section 3 along with our previous modeling results. Detailed two-dimensional microdefect simulations for point defects and point defect clusters are described in Sections 4 and 5 that demonstrate many of the features seen in crystal growth under moderate growth rates. In particular, the highlighted features are the formation of a vacancy-dominated core that leads to void formation, an annular ring near the crystal periphery with large self-interstitial super-saturation and eventual clustering, and the annular region that separates these regions, where the OSF-ring appears. The precise mechanism for the
formation of the OSF-ring is clear from the details of the vacancy profile within this region. These results give a context for understanding the proposed crystal growth of almost, defect free silicon under specific operating conditions.
2. Mathematical models of defect dynamics analysis The mathematical models for defect dynamics are developed in terms of microscopic transport theory by considering the equilibrium, diffusion, convection and reaction of defect species created by the combinations of vacancies, self-interstitials and impurities, such as oxygen. The development of these models is discussed here only for native point defects; the extension to include an impurity (boron) is described in [16] and the formulation for oxide precipitation is discussed in [17]. Microscopic balance equations are written for the vacancies (CV ðr; z; tÞ) and self-interstitials ðCI ðr; z; tÞÞ, where the position in the axisymmetric Czochralski crystal is described in cylindrical coordinates ðr; zÞ and time is given by t. The temperature field throughout the crystal is taken as Tðr; zÞ. The balance equations are given in [11,13] and include transport by Fickian diffusion, thermodiffusion, and crystal motion. Point defect loss and creation by recombination is model led by a law of mass action with the reaction rate given by kIV ðCI CV CIeq ðTÞCVeq ðTÞÞ;
ð1Þ Cieq ðTÞ
where kIV is the rate constant and is the equilibrium concentration of the ith species at temperature T. The loss of point defects to clusters is included as a source/sink term integrated over the entire cluster size distribution. The rate constant for recombination kIV ðTÞ is described by diffusion-limited reaction theory, coupled with a kinetic activation barrier; see [11,13]. The concentrations of small clusters are represented as discrete species of size n as Cn;X ðr; z; tÞ, where X represents either vacancy (V) or selfinterstitial (I) clusters. Large clusters of type X are represented by the continuous cluster size distribution written as fX ðr; z; tÞ. The discrete and continuous representations of the concentration are
100
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
matched at n ¼ nmatch , where the concentrations and fluxes are equated. Clusters are modeled as non-diffusing species. Discrete rate equations are written for the species Cn;X ðr; z; tÞ that account for convection caused by crystal motion and for growth and dissolution by monomer addition or subtraction from the cluster. The same physics is developed in the Fokker– Planck partial differential equation for the continuous portion of the size distribution fX ðr; z; tÞ. In both cases the rate constants for growth and dissolution rates are described by diffusion-limited theory with an activation barrier as 4pr2c DX DGn!nþ1;X kr;X ðn; r; z; tÞ ¼ exp ; ð2Þ q kT where rc is the radius of the cluster, DX is the diffusion coefficient of point defect X; q ¼ 0:235 nm is the lattice spacing, and k is the Boltzmann constant. The free energy change associated with incorporating monomer in the cluster is DGn!nþ1;X GfX ðn þ 1; TÞ GfX ðn; TÞ C1;X ðr; z; tÞ kT ln eq ; C1;X ðr; z; tÞ
ð3Þ
where G fX ðn; TÞ is the formation free energy of a cluster of size n of species X. In the calculations described here, the free energy of formation of the clusters is modeled using a very standard mesoscale phenomenological model for the growth or dissolution of spherical precipitates of type X with surface energy sX . This yields the functional form GfX ðn; TÞ ¼ aX sX n2=3 þ gX n þ G fX ð0Þ;
ð4Þ
where ðsX ; gX Þ are constants to be determined and aX is a constant that allows sX to be expressed as a surface energy as an equivalent spherical cluster. The term proportional to n takes into account the strain energy associated with the presence of the precipitate. The thermophysical properties for point defects used in the calculations presented here are the same as used in recent calculations [12,19,20], but differ slightly from those in [11,13]. The constant G fX ð0Þ is a correction to the phenomenological expression. For void formation, gV is set to zero because strain is not created in the lattice. We also
set gI ¼ 0; this value is consistent with parametric sensitivity analysis [12,20] for the nucleation temperature and with the failed attempts to image these defects by TEM. The surface energies are set to sV ¼ 1:0 J=m2 and sI ¼ 1:2 J=m2 . The models for self-interstitial formation do not include the evolution from spherical precipitates to stacking faults and, finally, perfect dislocation loops; including this detail will not alter appreciably the predictions of densities and sizes. The calculations presented here are based on numerical solution of the coupled set of partial differential equations for point defect dynamics, discrete rate equations and Fokker–Planck equations. The algorithm used is a newly developed combination [19] of the CC70 finite difference discretization of the n-space operator in the Fokker–Planck equation [21], Discontinuous Galerkin (DG) finite element discretizations in space and a semi-implicit time integration method. Because of the DG formulation [22] the discrete nonlinear algebraic equations that result are coupled only within each finite element and to neighboring elements through the convection term that models crystal motion; hence, the equations are readily solved in parallel using a very basic parallelization scheme. Parallel solution is essential because discretizations (typically 120 nodes in nspace, 60 elements in r and 280 elements in z) lead to over 20 million degrees-of-freedom in the finiteelement/finite-volume approximation.
3. The OSF-ring The OSF-ring has received an enormous amount of attention as a defect structure in CZgrown silicon. The ring is detectable by X-ray tomography of wafers after wet oxidation and is composed of extrinsic (interstitial-type) stacking fault loops lying in the h1 1 1i planes that grow anisotropically from small oxide precipitates [23]. It is thought that the self-interstitials are supplied by injection during the wet oxidation and collect on stacking faults that are already formed during oxide precipitation during crystal growth [24]. The structure of the OSF-ring is shown in the X-ray topograph shown in Fig. 2, which is taken from
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
101
Fig. 2. X-Ray topograph of radial cross-section of CZ crystal showing the vacancy-rich core (A), the OSF-ring (B), the nearly defect free-ring inside the OSF-ring (C), and the self-interstitial rich outer ring (D).
[24]. The ring separates the vacancy-rich core} which is populated by voids}from the selfinterstitial-rich outer annulus, which contains stacking faults. Interestingly, there is a finer structure to the microdefect distribution. Just outside the OSF-ring there is an annular region of almost microdefect-free crystal, as shown in Fig. 2. This region has received much attention [25,26] because its existence is an indication of the crystal growth conditions where almost perfect silicon can be grown. One of the goals of the modeling reported here is to reproduce this structure. The linkages between the radius of the OSF-ring and CZ processing conditions (pull rate and the temperature field in the crystal) have been well established; see the discussion in [13]. The most important empirical observation was reported by Dornberger et al. [14] who demonstrated that the radius of the ring could be correlated by the simple expression V ð5Þ ffi 1:34 103 cm2 =min K; GðROSF Þ where V is the pull rate and GðROSF Þ is the axial temperature gradient in the crystal measured at the melt/crystal interface and at the radial location r ¼ ROSF . Brown et al. [9] and then Sinno et al. [8] put forward the hypothesis that the approximate location of the OSF-ring could be predicted as the location of the neutral zone where neither point defect reached sufficient super-saturation to
allow nucleation and growth of the corresponding aggregates. For temperatures in the crystal much below the melting point (where the equilibrium point defect concentrations are much lower than the typical values in the crystal), this condition is expressed in terms of the excess concentration Dðr; zÞ CI ðr; zÞ CV ðr; zÞ as DðROSF ; zÞ ffi 0;
ð6Þ
where ROSF is the approximate location of the OSF-ring according to this criterion. Here it is assumed that the axial location is far enough from the solidification interface that the crystal has cooled sufficiently so that axial diffusion is no longer important. Sinno et al. [8] further assumed that the location r ¼ ROSF can be predicted solely from point defect dynamics, i.e. by solving the reaction–diffusion equations for point defects without accounting for cluster formation. This model of the separation between the vacancy-rich core and the self-interstitial-rich outer ring assumes that there is a spatial separation of point defect dynamics into distinct spatial regions in the crystal. This separation is shown schematically in Fig. 3. Sinno et al. [8] postulated there is a region near the melt/crystal interface where the high temperature of the crystal makes the point defects so mobile that the dynamics is dominated by diffusion and very rapid point defect recombination. As a result, intermediate, almost axially constant concentrations of point defects are established with the absolute
102
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 3. Schematic of regions for point defect reaction and formation of vacancy, self-interstitial and oxygen precipitates as a function of axial and radial location in the crystal.
concentrations depending on a delicate balance of axial diffusion and reaction at different radial locations across the crystal. As the crystal cools, these point defect concentrations may become super-saturated and nucleate clusters at nucleation temperatures that depend on the intermediate concentrations. As is experimentally well known, voids and self-interstitial defects nucleate at different temperatures. Actually, the nucleation temperatures are not constants, but depend on the intermediate point defect concentrations that supply the driving forces for clustering. The radial dependence of the temperature field and of these vacancy and self-interstitial concentrations cause radial variation of the nucleation region. Finally, at a yet lower temperature, oxide precipitates begin to form from the super-saturated oxygen concentration in the crystal. These precipitates preferentially form in a region of the crystal where the residual vacancy concentration in the wafer exists to supply the needed free volume
for precipitate growth according to the reaction Pn þ ð12 þ gI ÞSi þ O þ gV V $ Pnþ1 þ gI I þ stress energy; where Pn is a size of an n-mer oxide cluster and ðgI ; gV Þ are parameters computed to minimize the total energy of the oxide precipitate. Calculations presented here and in a series of papers have verified this picture, which is in qualitative agreement with the mechanisms proposed by Falster and Voronkov [25,26]. Using a combination of thermophysical properties predicted by atomistic simulations and by fitting entropic contributions to a single set of experimental data for the variation of the OSFring with crystal pull rate, Sinno et al. [8] were able to reproduce the expression in Eq. (5) using point defect simulations. These calculations also elucidated important features of the point defect structure. Most importantly, they demonstrated that point defect dynamics near the melt/ crystal interface are dominated by very rapid
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
recombination that is quenched at even small axial distances (z5R) because of crystal cooling. For this limit, an appropriate dimensionless form of the equations leads to boundary-layer analysis in which the variation of the vacancy and selfinterstitial concentrations in this layer are predicted solely as a function of a linear approximation to the axial temperature gradient GðrÞ and the thermophysical properties of the point defects. Using these concentration profiles in the condition for the neutral ring, Eq. (5), yielded the closedform expression: V GðRD¼0 Þ DI ðTm ÞCIeq ðTm ÞEI* DV ðTm ÞCVeq ðTm ÞEV* ¼ Tm ðCVeq ðTm Þ CIeq ðTm ÞÞ ¼ K0 103 cm2 min1 K1 ;
ð7Þ
where EX* are constants computed from the leading order asymptotic solution (see [8]) and DX ðTm Þ and CXeq ðTm Þ are the point defect diffusivities and equilibrium concentrations evaluated at the melting temperature Tm . Using the thermo-
103
physical properties in [8], Eq. (7) gives K0 ¼ 1:38 103 cm2 min1 K 1 ; using the properties used in this paper and in [12] gives K0 ¼ 1:43 103 cm2 min1 K 1 . Although the point defect simulations and analysis that were first reported in [8] gave a quantitative description for the evolution of OSFring with processing conditions, several features were not explained. These included the reason for the preferential formation of oxide precipitates near the region where DðRD¼0 ; zÞ ffi 0, the finite width of the ring and the appearance of the nearly defect free region outside the OSF-ring. These issues are clear in the two-dimensional microdefect calculations reported below. 4. Two-dimensional simulations of defect dynamics We report calculations for a prototype CZ-crystal with radius R ¼ 7:8 cm grown in a furnace that yields a specific temperature profile in the crystal. The temperature field is shown in Fig. 4
Fig. 4. Temperature field used in defect dynamics simulations; the axial temperature gradient G at the interface also is plotted for the temperature field.
104
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
and has variation in both the axial and radial dimensions. In this and other figures, the coordinates ðr; zÞ have been made dimensionless with the length scale L ¼ 10 cm, so that the radius of the crystal corresponds to r ¼ 0:78. The calculations reported here do not include oxygen dynamics or oxide formation. We mimic the first stages of this process by stopping vacancy precipitation into voids at temperatures below T ¼ 9508C, where we believe these vacancies would preferentially react with oxygen to begin the precipitation process. Prediction of point defect dynamics near the surface of the crystal depends sensitively on
the description of the interactions of point defects on the surface. We use the model described by Mori and Brown [12,20] where vacancies and self-interstitials react by surface recombination according to a mixed or Robin boundary condition. A sample steady-state solution is shown in Fig. 5 computed for the pull rate V ¼ 0:6 mm=min. Several features of the results are obvious. These include the rapid decrease in the vacancy and self-interstitial concentrations near the melt/crystal interface (A), the intermediate high concentration of vacancies (approximately
Fig. 5. Simulation results for steady-state crystal growth at V ¼ 0:6 mm=min showing (a) vacancy concentration CV ðr; zÞ, (b) selfinterstitial density CI ðr; zÞ, (c) total void cluster density NV ðr; zÞ (>50 nm in diameter), and self-interstitial cluster density NI ðr; zÞ(>50 nm in diameter). Size distributions for voids and self-interstitial clusters at the top of the crystal are shown in an inset.
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
4 1013 cm3 ) near the center of the crystal and the high self-interstitial concentrations near the periphery (B); these two features also are shown in Fig. 7a below. The location DðRD¼0 ; zÞ ffi 0 occurs at approximately r ¼ 0:4, as shown at (C). At the top of the crystal, both the vacancy and selfinterstitial concentrations have decreased because of precipitation, as shown by the explosive growth in the total cluster concentrations NV ðr; zÞ and NI ðr; zÞ (defined as the sum of all clusters of size greater that 50 nm). In these simulations, the predicted onset value for vacancy nucleation is approximately 10708C and is at the lower end of the range, 10708C4T411508C, seen in experiments [27–29]. The final cluster density predicted along the centerline is approximately NV ¼ 4 105 cm3 and is well within the range seen in experiments. The temperature for self-interstitial nucleation is predicted to be 9408C, as compared to the experimental range 9508C4T410308C seen in experiments [30,31]. The size distributions for
105
these defects predicted at the top of the crystal are shown in the inset in Fig. 5. The variations of the point defect profiles with pull rate V are demonstrated in Figs. 6 and 7. Decreasing the pull rate causes the region of void formation to shrink toward the centerline of the crystal and for V ¼ 0:5 mm/min this region has disappeared and no voids are nucleated. Correspondingly, self-interstitial super-saturation occurs in a larger portion of the crystal and more selfinterstitial clusters are produced. Alternatively, increasing pull rate extends the region of void formation to a larger radial core of the crystal and pushes the region of self-interstitials toward the periphery. For V ¼ 0:9 mm/ min, the self-interstitial rich region has been moved so close to the surface that diffusion to and recombination with at the surface of point defects has severely decreased the maximum in the intermediate CI concentration and self-interstitial clusters are not observed.
Fig. 6. Vacancy fields predicted at the pull rate values of 0.5 mm/min, 0.7 mm/min, and 0.9 mm/min.
106
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 8. Variation of V=GðRD¼0 Þ and V=GðRCV;max Þ with pull rate. The asymptotic result from Eq. (7) is shown for comparison.
This appearance of the peak suggests another definition for the location of the middle of the OSF-ring, as the location of the maximum in the residual vacancy profile after void nucleation; we define this value as RCV;max . This leads to ratio V ð8Þ ¼ K0max cm2 min1 K 1 GðRCV;max Þ
Fig. 7. Radial profiles of CV and CI computed (a) at z ¼ 2:8 cm (intermediate concentrations) and (b) near the top of the crystal (residual vacancy concentration) for several pull rates. The location of the OSF-ring corresponds to the region with high residual vacancy concentration in (b).
5. Residual vacancy concentrations, the osf-ring and almost perfect silicon An important feature of Fig. 7b is the appearance of the peak in the residual vacancy concentration at a particular value of the radius, which increases with pull rate.
as the more precise interpretation of the experiments of Dornberger et al. [14]. The expression Eq. (8) is plotted from our simulations on Fig. 8 and gives 1:50 103 4K0max 41:55 103 , compared to the value of K0 ¼ 1:43 103 computed from the asymptotic expression Eq. (7) and the range 1:28 103 4K0 41:35 103 extracted by plotting the value of V=GðRD¼0 Þ computed from the intermediate concentrations. The variation of K0 with both V=GðRD¼0 Þ and C=GðRCV;max Þ predicted from the simulations is partially a result of the realistic modeling of two-dimensional effects at the crystal surface introduced in the simulations. As shown in Fig. 6 and, certainly in the radial concentration profiles in Fig. 7a, the self-interstitial concentrations are affected by the presence of the crystal surface for pull rates as low as 0.7 mm/ min. Only in crystal growth of larger diameter crystals (800 diameter and above) would one expect
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
an appreciable range of V/G where K0 may be nearly constant. The values of K0 and K max are well within the 0 range of error inherent in estimating the interfacial temperature gradients G in the experiments. Moreover, the fact that the maximum in the residual vacancy concentration moves approximately with pull rate according to V=G suggests that this location can be recovered from an asymptotic analysis as well. However, such an analysis must include both the formation of voids and the effect of variations in the intermediate vacancy concentration. Analysis of the two-dimensional simulations also shows the fine-scale structure of the nearly defect-free region just outside of the OSF-ring. This structure is best seen by the plots in Fig. 9 of the radial structure of the total defects concentrations NV ðr; zÞ and NI ðr; zÞ computed near the top of the crystal for several pull rates. Also shown on these figures are the locations r ¼ RCV;max ðVÞ and r ¼ RD¼0 ðVÞ taken from the simulation for
107
V ¼ 0:6 mm/min. The region of low defect concentration is clearly visible near r ¼ Rfree ðVÞ and RCV;max 4Rfree 4RD¼0 . As with other markers for the radial structure, Rfree ðVÞvaries with V and according to V=G ðRfree Þ ¼ K0free cm2 =minK, where K0free is approximately equal to K0 . Again this result suggests that an asymptotic analysis can predict this dependence. The time-dependent microdefect dynamics simulations are useful tools to predict the spatial and temporal evolution of the microdefect structure with changes in operating conditions, such as the crystal pull rate V. The predicted evolution of the time-dependent dynamics of point defects and their clusters are shown in Fig. 10 for changing pull rate, but with a constant temperature field. The OSF-ring and IV boundaries separate as they approach the crystal surface with increasing pull rate. The number density of interstitial clusters decreases well as the IV neutral boundary comes close to the crystal surface and self-interstitials leave the crystal by diffusion and surface recombination. A lifetime map also is shown for (an experimental measure of the residual vacancy concentration in silicon) a horizontally sectioned CZ crystal in which a qualitatively similar growth rate schedule was used. The qualitative agreement in the two-dimensional structure of the boundary separating the void and self-interstitial regions is very good. More detailed comparison will have to await precise experiments and matching simulations.
6. Summary
Fig. 9. Radial profiles of total defect concentrations ðNV ðr; zÞ; NI ðr; zÞÞ computed at an axial position near the top of the crystal for pull rates (V) of 0.6 mm/min, 0.7 mm/min and 0.9 mm/min. Most significant is the prediction of an annular region of almost microdefect defect-free (of size greater than 50 nm) silicon as shown by the point (P) for V ¼ 0:7 mm/min, the locations of ðRD¼0 ; RCV;max ; Rfree Þ also are indicated for this pull rate.
The simulations reported here and in other references demonstrate the impact of quantitative analysis of defect dynamics on understanding microdefect formation during silicon crystal growth and subsequent processing. An enormous number of the patterns of microdefects seen during silicon crystal growth are semi-quantitatively predicted using self-consistent models of the equilibrium, transport, reaction and clustering of point defects. Quantitative agreement seems achievable, but is a strong function of accuracy of the thermophysical properties for point defects
108
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 10. Simulation results for step change in crystal growth rate showing (a) vacancy density, and (b) lifetime map for crystal grown with similar step change in pull rate. Life time map of silicon crystal grown under similar conditions.
and of the level of detail incorporated into the simulations. The simulations reported here only described the interactions of point defects. Incorporation to include impurities, such as boron, is feasible [16]. More importantly, the extension of the analysis to include oxide precipitate formation has been accomplished [17]. When one-dimensional simulations including oxide precipitation are performed for temperature fields corresponding to the radial location of the maximum in residual vacancy concentration, r ¼ RCv;max , a measurable increase in oxygen precipitation is predicted during crystal cooling; these oxide precipitates are postulated to be the precursors of the defects in the OSF-ring. More interestingly, these defects are easily dissolved during a high temperature wafer annealing, leaving an almost defect-free (low values of NV and NI ) region in the crystal. Capitalizing on the appearance of the defect-free crystal region is difficult. To grow crystals with this
property across a wide portion of the radius would require designing a CZ hot-zone with an approximately radially uniform crystal temperature tailored so that defect clustering is avoided, independent of changes in pull rate and melt volume. This is a difficult challenge and it points to the need to couple the microdefect simulations to accurate macroscale global simulations of heat transfer and melt convection in commercial-scale CZ systems. The development of such simulations is underway [32]. The coupling through the temperature field in the crystal of these simulations with defect dynamics simulations will make the vision of simulation-aided defect engineering of large-scale CZ systems a reality [33].
Acknowledgements We thank MEMC Corporation, Shin-Etsu Hanotai Company and Wacker Siltronic for
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
financial support of this research. We are indebted to Erich Dornberger and Christian Hoess for the lifetime map of the crystal shown in Fig. 10.
References [1] A.J.R. de Kock, J. Wijgert, J. Crystal Growth 49 (1980) 718. [2] H. Yamagishi, I. Fusegawa, N. Fujimaki, M. Katayama, Semicond. Sci. Technol. 7 (1992) A135. [3] J. Ryuta, E. Morita, T. Tanaka, Y. Shimanuki, Jpn. J. Appl. Phys. 29 (1990) 1947. [4] M. Hourai, T. Nagashima, E. Kajita, S. Miki, T. Shigematsu, M. Okiu, J. Electrochem. Soc. 142 (1995) 3193. [5] T.S. Plaskett, Trans. Met. Soc. AIME, 233 (1965) 809. [6] H. Fo¨ll, U. Go¨sele, B.O. Kolbesen, J. Crystal Growth 52 (1981) 907. [7] T. Abe, H. Harada, MRS Symp. Proc. 14 (1983) 1. [8] T. Sinno, R.A. Brown, W. von Ammon, E. Dornberger, J. Electrochem. Soc. 145 (1998) 302. [9] V.V. Voronkov, J. Crystal Growth 59 (1982) 625. [10] R.A. Brown, D. Maroudas, T. Sinno, J. Crystal Growth 137 (1994) 12. [11] T. Sinno, R.A. Brown, J. Electrochem. Soc. 146 (1999) 2300. [12] T. Mori, Doctoral dissertation, Massachusetts Institute of Technology, Cambridge, MA, 2000. [13] T. Sinno, E. Dornberger, W. von Ammon, R.A. Brown, F. Dupret, Mater. Sci. Eng. 243 (2000) 1. [14] E. Dornberger, D. Gra¨f, M. Suhren, U. Lambert, P. Wagner, F. Dupret, W. von Ammon, J. Crystal Growth 180 (1997) 343. [15] T. Ueki, M. Itsumi, T. Takeda, Jpn. J. Appl. Phys. 38 (1999) 5695. [16] T. Sinno, H. Susanto, R.A. Brown, Appl. Phys. Lett. 75 (1999) 1544.
109
[17] Z. Wang, R.A. Brown, J. Crys. Growth (2000) submitted. [18] A.K. Kapila, Asymptotic Treatment of Chemically Reacting Systems, Pitman, London, 1983. [19] T. Mori, R.A. Brown, Chem. Eng. Sci. (2001) submitted. [20] T. Mori, Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, 2000. [21] J.S. Chang, G. Cooper, J. Comp. Phys. 6 (1970) 1. [22] B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods, Springer, New York, 1999. [23] K. Sueoka, M. Akatsuk, K. Nishihara, T. Yamamoto, S. Kobayashi, Mater. Sci. Forum 196–201 (1995) 1737. [24] M. Hasebe, Y. Takeoka, A. Shinoyama, S. Naito, J. Appl. Phys. 28 (1989) 1999. [25] V.V. Voronkov, R. Falster, J. Crystal Growth 194 (1998) 76. [26] V.V. Voronkov, R. Falster, J. Crystal Growth 204 (1999) 462. [27] M. Hourai, M. Nagashima, T. Kajita, E. Miki, S. Shigematsu, T. Okui, J. Electrochem Soc. 142 (1995) 3193. [28] K. Takano, K. Kitagawa, E. Iino, M. Kimura, H. Yamagishi, Mater. Sci. Forum 196 (1995) 1707. [29] R. Falster, V.V. Voronkov, J.C. Holzer, S. Markgraf, S. Macquaid, L. Mule’Stagno, In: H.R. Huff, H. Tsuya, U. Go¨sele, (Eds.), Semiconductor Silicon 1998, Electrochemical Society, Pennington, 1998. [30] M. Okiu, T. Tananka, T. Kanada, T. Ono, Ouyou Butsiri 66 (1997) 707. [31] T Saishouji, L. Nakamura, H. Nakajima, T. Yokoyama, F. Ishikawa, J. Tomioka, In: C.L. Claeys, P. RaiChoudhury, M. Watanabe, P. Stallhofer, H.J. Dawson, (Eds.), High Purity Silicon V, Electrochemical Society, Pennington, 1998. [32] A. Lipchin, R.A. Brown, J. Crystal Growth 216 (2000) 192–203. [33] R.A. Brown, Z. Wang, Y.C. Won, and A. Lipchin, J. Crystal Growth (2001) Invited paper, submitted.
Engineering analysis of microdefect formation during silicon crystal growth Robert A. Brown*, Zhihong Wang, Tatsuo Mori Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave 3-208, Cambridge, MA 02139-4307, USA
Abstract The quality of single crystalline silicon wafers used as substrates for microelectronic devices is measured in terms of the type, size and number density of microdefects formed during crystal growth and subsequent processing. Native point defects}vacancies and self-interstitials}and impurities, such as oxygen and carbon, in the silicon crystal convect, diffuse, react and aggregate, driven by species super-saturation, to form defect structures at densities that are microscopically visible. We model defect dynamics using diffusion–reaction theory in which defect types are represented as chemical species and have concentrations that are governed by conservation laws, chemical thermodynamics and kinetics, using self-consistent values of equilibrium, transport and kinetic parameters for point defects and impurities. We refer to this approach as crystal defect dynamics analysis. The purpose of this paper is to demonstrate that this approach leads to a wide range of predictions that are consistent with experiments.The results presented here focus on the two-dimensional predictions of microdefect structures under growth conditions where both voids (cluster of vacancies) and self-interstitial aggregates are seen in regions of the crystal separated by the oxidation-induced stacking fault ring (OSF-ring). The location and structure of the OSF-ring is predicted in terms of point defect dynamics near the melt/crystal interface, which results in an annular region with almost balanced point defect concentrations (centered at r RD¼0 ), and in terms of a peak (centered at r RCV;max ) in the residual vacancy concentration after cluster formation. This residual vacancy concentration facilitates oxide precipitation during crystal growth, which seeds stacking fault formation during subsequent oxidation. Calculations show that RCV;max 4RD¼0 with both values scaling with V=GðRÞ, where V is the crystal pull rate and GðRÞ is the axial temperature gradient at r ¼ R along the melt/crystal interface. An annulus just outside the OSF-ring (centered at r Rfree ; RCV;max 4Rfree 4RD¼0 ) is identified where almost no microdefects are present, thereby identifying the temperature field that leads to almost perfect silicon. # 2001 Elsevier Science B.V. All rights reserved. Keywords: A1. Computer simulation; A1. Defects; B2. Semiconducting silicon
1. Introduction In the last decade there has been tremendous progress in the understanding of the cause for the *Corresponding author. Tel.: +1-617-253-5726; fax: +1617-253-8549. E-mail address: [email protected] (R.A. Brown).
myriad of crystallographic imperfections that are present in the otherwise perfect single-crystal silicon substrates grown by the Czochralski (CZ) crystal growth method. Today it is becoming well understood that the microdefects observed in asgrown silicon boules and the precipitation of oxygen and other impurities during growth and wafer processing is controlled by the dynamics of
0022-0248/01/$ - see front matter # 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 1 ) 0 0 8 2 5 - 9
98
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
intrinsic point defects}vacancies and self-interstitials}which diffuse, react and aggregate to form the structures observed by etching methods [1,2], particle counting [3], laser scattering [4] and other techniques. A new approach has been developed with the aim of connecting the dynamics of point defects and microdefect evolution with processing conditions during crystal growth and subsequent wafer processing; we call this approach crystal defect dynamics analysis. The role of intrinsic point defects in microdefect formation has long been postulated in the formation of microdefects such as the self-interstitial aggregates known as A- and B-defects [5,6] and the vacancy dominated D-defects in float zone silicon [7]. The transition between vacancy-rich defects at the core of a Czochralski crystal and selfinterstitial dominated defects near the periphery is shown in the X-ray topograph of the axial section of a CZ crystal in Fig. 1, taken from [8]. The experiment shows the change in the defect structure with variation in pull rate. The sharp band separating voids and stacking faults is denoted and responds to pull rate changes, moving inward for decreasing crystal pull rate and outward with increasing pull rate.
Fig. 1. X-ray topograph of crystallographic defects in axial cross-section of CZ crystal. The crystal was grown with variations in the pull rate, as shown by the scale on the left, The migration of the OSF-ring with pull rate is clearly seen.
Voronkov [9] was the first to connect theoretically point defect dynamics and crystal growth operating conditions, such as crystal pull rate V and temperature field, with observable transitions in microdefect structures during crystal growth. His approach has pioneered comprehensive modeling of point defect dynamics [8,10–13], which is the subject of this paper. The analysis used in this approach is built on the hypothesis that microdefects in silicon crystals are the result of the transport, equilibrium and reaction of intrinsic point defects and impurities. The reactions may be to either annihilate defects, as in the case of selfinterstitial and vacancy recombination, or to nucleate and aggregate defects to form larger, more complex species. Large clusters form into microscopically observable defects driven by large point defect (and impurity) super-saturations that result as the crystal cools and during wafer processing. Subtle changes in the temperature field during processing [8] and changes in point defect dynamics caused by interactions with impurities, such as boron and carbon [14–16] swing the balance from vacancies to self-interstitials and change the final microdefect structure. The success of point defect dynamics for the prediction of defect structures in silicon has been phenomenal. The predictions include understanding the sensitivity to the crystal temperature field and pull rate of important features of the microdefects in the crystal. These features include; the boundary between the vacancy- and selfinterstitial-rich regions in the crystal, the nucleation temperature, the defect number density and average observable size of microvoids in vacancydominated crystal growth and the onset temperature and density of self-interstitial defects. Most recently, the predictions have been extended to many of the features of oxygen precipitation observed during crystal growth and wafer annealing [17]. Most remarkably, a self-consistent set of predictions that are in reasonable agreement with experiment can be obtained using thermophysical properties for point defect equilibrium and transport grounded in atomistic simulations and fit to experimental data for the dynamics of the oxidation-induced stacking fault ring (OSF-ring), as described elsewhere [8,13].
99
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
The defect dynamics models used to describe these processes are built on simple concepts for transport and diffusion-controlled reaction for dilute species in a crystalline solid. The formulation is made more complicated by the need to describe clusters of point defects and impurities that range in size from dimers, trimers, etc. to clusters as large as 1400 nm in diameter, which contain Oð1010 Þ point defects! This range is achieved by coupling the discrete rate equations that describe the reactions of small clusters with the Fokker–Planck or continuous rate description for the dynamics of large species [11]. The formulation of the defect dynamics models is described in detail elsewhere and is only reviewed in Section 2 [11–13]. These formulations lead to extremely complex combinations of diffusion– reaction equations for the defect species representing monomeric species and aggregates. These equations are similar to equations describing classical chemical reactions [18]. We have developed novel numerical methods for the accurate and efficient parallel solution of these equations for both steady and transient dynamics [19]. Solutions and experimental evidence have suggested asymptotic analysis for specific situations that has provided valuable insight into scaling [8]. The purpose of this paper is to demonstrate the application of defect dynamics analysis for several specific problems that arise in Czochralski crystal growth. Of specific interest are the details of the formation of the OSF-ring in large-scale crystals. We use the analysis of this microdefect structure as the basis for the presentation that follows. The experimental features of the OSF-ring are reviewed in Section 3 along with our previous modeling results. Detailed two-dimensional microdefect simulations for point defects and point defect clusters are described in Sections 4 and 5 that demonstrate many of the features seen in crystal growth under moderate growth rates. In particular, the highlighted features are the formation of a vacancy-dominated core that leads to void formation, an annular ring near the crystal periphery with large self-interstitial super-saturation and eventual clustering, and the annular region that separates these regions, where the OSF-ring appears. The precise mechanism for the
formation of the OSF-ring is clear from the details of the vacancy profile within this region. These results give a context for understanding the proposed crystal growth of almost, defect free silicon under specific operating conditions.
2. Mathematical models of defect dynamics analysis The mathematical models for defect dynamics are developed in terms of microscopic transport theory by considering the equilibrium, diffusion, convection and reaction of defect species created by the combinations of vacancies, self-interstitials and impurities, such as oxygen. The development of these models is discussed here only for native point defects; the extension to include an impurity (boron) is described in [16] and the formulation for oxide precipitation is discussed in [17]. Microscopic balance equations are written for the vacancies (CV ðr; z; tÞ) and self-interstitials ðCI ðr; z; tÞÞ, where the position in the axisymmetric Czochralski crystal is described in cylindrical coordinates ðr; zÞ and time is given by t. The temperature field throughout the crystal is taken as Tðr; zÞ. The balance equations are given in [11,13] and include transport by Fickian diffusion, thermodiffusion, and crystal motion. Point defect loss and creation by recombination is model led by a law of mass action with the reaction rate given by kIV ðCI CV CIeq ðTÞCVeq ðTÞÞ;
ð1Þ Cieq ðTÞ
where kIV is the rate constant and is the equilibrium concentration of the ith species at temperature T. The loss of point defects to clusters is included as a source/sink term integrated over the entire cluster size distribution. The rate constant for recombination kIV ðTÞ is described by diffusion-limited reaction theory, coupled with a kinetic activation barrier; see [11,13]. The concentrations of small clusters are represented as discrete species of size n as Cn;X ðr; z; tÞ, where X represents either vacancy (V) or selfinterstitial (I) clusters. Large clusters of type X are represented by the continuous cluster size distribution written as fX ðr; z; tÞ. The discrete and continuous representations of the concentration are
100
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
matched at n ¼ nmatch , where the concentrations and fluxes are equated. Clusters are modeled as non-diffusing species. Discrete rate equations are written for the species Cn;X ðr; z; tÞ that account for convection caused by crystal motion and for growth and dissolution by monomer addition or subtraction from the cluster. The same physics is developed in the Fokker– Planck partial differential equation for the continuous portion of the size distribution fX ðr; z; tÞ. In both cases the rate constants for growth and dissolution rates are described by diffusion-limited theory with an activation barrier as 4pr2c DX DGn!nþ1;X kr;X ðn; r; z; tÞ ¼ exp ; ð2Þ q kT where rc is the radius of the cluster, DX is the diffusion coefficient of point defect X; q ¼ 0:235 nm is the lattice spacing, and k is the Boltzmann constant. The free energy change associated with incorporating monomer in the cluster is DGn!nþ1;X GfX ðn þ 1; TÞ GfX ðn; TÞ C1;X ðr; z; tÞ kT ln eq ; C1;X ðr; z; tÞ
ð3Þ
where G fX ðn; TÞ is the formation free energy of a cluster of size n of species X. In the calculations described here, the free energy of formation of the clusters is modeled using a very standard mesoscale phenomenological model for the growth or dissolution of spherical precipitates of type X with surface energy sX . This yields the functional form GfX ðn; TÞ ¼ aX sX n2=3 þ gX n þ G fX ð0Þ;
ð4Þ
where ðsX ; gX Þ are constants to be determined and aX is a constant that allows sX to be expressed as a surface energy as an equivalent spherical cluster. The term proportional to n takes into account the strain energy associated with the presence of the precipitate. The thermophysical properties for point defects used in the calculations presented here are the same as used in recent calculations [12,19,20], but differ slightly from those in [11,13]. The constant G fX ð0Þ is a correction to the phenomenological expression. For void formation, gV is set to zero because strain is not created in the lattice. We also
set gI ¼ 0; this value is consistent with parametric sensitivity analysis [12,20] for the nucleation temperature and with the failed attempts to image these defects by TEM. The surface energies are set to sV ¼ 1:0 J=m2 and sI ¼ 1:2 J=m2 . The models for self-interstitial formation do not include the evolution from spherical precipitates to stacking faults and, finally, perfect dislocation loops; including this detail will not alter appreciably the predictions of densities and sizes. The calculations presented here are based on numerical solution of the coupled set of partial differential equations for point defect dynamics, discrete rate equations and Fokker–Planck equations. The algorithm used is a newly developed combination [19] of the CC70 finite difference discretization of the n-space operator in the Fokker–Planck equation [21], Discontinuous Galerkin (DG) finite element discretizations in space and a semi-implicit time integration method. Because of the DG formulation [22] the discrete nonlinear algebraic equations that result are coupled only within each finite element and to neighboring elements through the convection term that models crystal motion; hence, the equations are readily solved in parallel using a very basic parallelization scheme. Parallel solution is essential because discretizations (typically 120 nodes in nspace, 60 elements in r and 280 elements in z) lead to over 20 million degrees-of-freedom in the finiteelement/finite-volume approximation.
3. The OSF-ring The OSF-ring has received an enormous amount of attention as a defect structure in CZgrown silicon. The ring is detectable by X-ray tomography of wafers after wet oxidation and is composed of extrinsic (interstitial-type) stacking fault loops lying in the h1 1 1i planes that grow anisotropically from small oxide precipitates [23]. It is thought that the self-interstitials are supplied by injection during the wet oxidation and collect on stacking faults that are already formed during oxide precipitation during crystal growth [24]. The structure of the OSF-ring is shown in the X-ray topograph shown in Fig. 2, which is taken from
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
101
Fig. 2. X-Ray topograph of radial cross-section of CZ crystal showing the vacancy-rich core (A), the OSF-ring (B), the nearly defect free-ring inside the OSF-ring (C), and the self-interstitial rich outer ring (D).
[24]. The ring separates the vacancy-rich core} which is populated by voids}from the selfinterstitial-rich outer annulus, which contains stacking faults. Interestingly, there is a finer structure to the microdefect distribution. Just outside the OSF-ring there is an annular region of almost microdefect-free crystal, as shown in Fig. 2. This region has received much attention [25,26] because its existence is an indication of the crystal growth conditions where almost perfect silicon can be grown. One of the goals of the modeling reported here is to reproduce this structure. The linkages between the radius of the OSF-ring and CZ processing conditions (pull rate and the temperature field in the crystal) have been well established; see the discussion in [13]. The most important empirical observation was reported by Dornberger et al. [14] who demonstrated that the radius of the ring could be correlated by the simple expression V ð5Þ ffi 1:34 103 cm2 =min K; GðROSF Þ where V is the pull rate and GðROSF Þ is the axial temperature gradient in the crystal measured at the melt/crystal interface and at the radial location r ¼ ROSF . Brown et al. [9] and then Sinno et al. [8] put forward the hypothesis that the approximate location of the OSF-ring could be predicted as the location of the neutral zone where neither point defect reached sufficient super-saturation to
allow nucleation and growth of the corresponding aggregates. For temperatures in the crystal much below the melting point (where the equilibrium point defect concentrations are much lower than the typical values in the crystal), this condition is expressed in terms of the excess concentration Dðr; zÞ CI ðr; zÞ CV ðr; zÞ as DðROSF ; zÞ ffi 0;
ð6Þ
where ROSF is the approximate location of the OSF-ring according to this criterion. Here it is assumed that the axial location is far enough from the solidification interface that the crystal has cooled sufficiently so that axial diffusion is no longer important. Sinno et al. [8] further assumed that the location r ¼ ROSF can be predicted solely from point defect dynamics, i.e. by solving the reaction–diffusion equations for point defects without accounting for cluster formation. This model of the separation between the vacancy-rich core and the self-interstitial-rich outer ring assumes that there is a spatial separation of point defect dynamics into distinct spatial regions in the crystal. This separation is shown schematically in Fig. 3. Sinno et al. [8] postulated there is a region near the melt/crystal interface where the high temperature of the crystal makes the point defects so mobile that the dynamics is dominated by diffusion and very rapid point defect recombination. As a result, intermediate, almost axially constant concentrations of point defects are established with the absolute
102
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 3. Schematic of regions for point defect reaction and formation of vacancy, self-interstitial and oxygen precipitates as a function of axial and radial location in the crystal.
concentrations depending on a delicate balance of axial diffusion and reaction at different radial locations across the crystal. As the crystal cools, these point defect concentrations may become super-saturated and nucleate clusters at nucleation temperatures that depend on the intermediate concentrations. As is experimentally well known, voids and self-interstitial defects nucleate at different temperatures. Actually, the nucleation temperatures are not constants, but depend on the intermediate point defect concentrations that supply the driving forces for clustering. The radial dependence of the temperature field and of these vacancy and self-interstitial concentrations cause radial variation of the nucleation region. Finally, at a yet lower temperature, oxide precipitates begin to form from the super-saturated oxygen concentration in the crystal. These precipitates preferentially form in a region of the crystal where the residual vacancy concentration in the wafer exists to supply the needed free volume
for precipitate growth according to the reaction Pn þ ð12 þ gI ÞSi þ O þ gV V $ Pnþ1 þ gI I þ stress energy; where Pn is a size of an n-mer oxide cluster and ðgI ; gV Þ are parameters computed to minimize the total energy of the oxide precipitate. Calculations presented here and in a series of papers have verified this picture, which is in qualitative agreement with the mechanisms proposed by Falster and Voronkov [25,26]. Using a combination of thermophysical properties predicted by atomistic simulations and by fitting entropic contributions to a single set of experimental data for the variation of the OSFring with crystal pull rate, Sinno et al. [8] were able to reproduce the expression in Eq. (5) using point defect simulations. These calculations also elucidated important features of the point defect structure. Most importantly, they demonstrated that point defect dynamics near the melt/ crystal interface are dominated by very rapid
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
recombination that is quenched at even small axial distances (z5R) because of crystal cooling. For this limit, an appropriate dimensionless form of the equations leads to boundary-layer analysis in which the variation of the vacancy and selfinterstitial concentrations in this layer are predicted solely as a function of a linear approximation to the axial temperature gradient GðrÞ and the thermophysical properties of the point defects. Using these concentration profiles in the condition for the neutral ring, Eq. (5), yielded the closedform expression: V GðRD¼0 Þ DI ðTm ÞCIeq ðTm ÞEI* DV ðTm ÞCVeq ðTm ÞEV* ¼ Tm ðCVeq ðTm Þ CIeq ðTm ÞÞ ¼ K0 103 cm2 min1 K1 ;
ð7Þ
where EX* are constants computed from the leading order asymptotic solution (see [8]) and DX ðTm Þ and CXeq ðTm Þ are the point defect diffusivities and equilibrium concentrations evaluated at the melting temperature Tm . Using the thermo-
103
physical properties in [8], Eq. (7) gives K0 ¼ 1:38 103 cm2 min1 K 1 ; using the properties used in this paper and in [12] gives K0 ¼ 1:43 103 cm2 min1 K 1 . Although the point defect simulations and analysis that were first reported in [8] gave a quantitative description for the evolution of OSFring with processing conditions, several features were not explained. These included the reason for the preferential formation of oxide precipitates near the region where DðRD¼0 ; zÞ ffi 0, the finite width of the ring and the appearance of the nearly defect free region outside the OSF-ring. These issues are clear in the two-dimensional microdefect calculations reported below. 4. Two-dimensional simulations of defect dynamics We report calculations for a prototype CZ-crystal with radius R ¼ 7:8 cm grown in a furnace that yields a specific temperature profile in the crystal. The temperature field is shown in Fig. 4
Fig. 4. Temperature field used in defect dynamics simulations; the axial temperature gradient G at the interface also is plotted for the temperature field.
104
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
and has variation in both the axial and radial dimensions. In this and other figures, the coordinates ðr; zÞ have been made dimensionless with the length scale L ¼ 10 cm, so that the radius of the crystal corresponds to r ¼ 0:78. The calculations reported here do not include oxygen dynamics or oxide formation. We mimic the first stages of this process by stopping vacancy precipitation into voids at temperatures below T ¼ 9508C, where we believe these vacancies would preferentially react with oxygen to begin the precipitation process. Prediction of point defect dynamics near the surface of the crystal depends sensitively on
the description of the interactions of point defects on the surface. We use the model described by Mori and Brown [12,20] where vacancies and self-interstitials react by surface recombination according to a mixed or Robin boundary condition. A sample steady-state solution is shown in Fig. 5 computed for the pull rate V ¼ 0:6 mm=min. Several features of the results are obvious. These include the rapid decrease in the vacancy and self-interstitial concentrations near the melt/crystal interface (A), the intermediate high concentration of vacancies (approximately
Fig. 5. Simulation results for steady-state crystal growth at V ¼ 0:6 mm=min showing (a) vacancy concentration CV ðr; zÞ, (b) selfinterstitial density CI ðr; zÞ, (c) total void cluster density NV ðr; zÞ (>50 nm in diameter), and self-interstitial cluster density NI ðr; zÞ(>50 nm in diameter). Size distributions for voids and self-interstitial clusters at the top of the crystal are shown in an inset.
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
4 1013 cm3 ) near the center of the crystal and the high self-interstitial concentrations near the periphery (B); these two features also are shown in Fig. 7a below. The location DðRD¼0 ; zÞ ffi 0 occurs at approximately r ¼ 0:4, as shown at (C). At the top of the crystal, both the vacancy and selfinterstitial concentrations have decreased because of precipitation, as shown by the explosive growth in the total cluster concentrations NV ðr; zÞ and NI ðr; zÞ (defined as the sum of all clusters of size greater that 50 nm). In these simulations, the predicted onset value for vacancy nucleation is approximately 10708C and is at the lower end of the range, 10708C4T411508C, seen in experiments [27–29]. The final cluster density predicted along the centerline is approximately NV ¼ 4 105 cm3 and is well within the range seen in experiments. The temperature for self-interstitial nucleation is predicted to be 9408C, as compared to the experimental range 9508C4T410308C seen in experiments [30,31]. The size distributions for
105
these defects predicted at the top of the crystal are shown in the inset in Fig. 5. The variations of the point defect profiles with pull rate V are demonstrated in Figs. 6 and 7. Decreasing the pull rate causes the region of void formation to shrink toward the centerline of the crystal and for V ¼ 0:5 mm/min this region has disappeared and no voids are nucleated. Correspondingly, self-interstitial super-saturation occurs in a larger portion of the crystal and more selfinterstitial clusters are produced. Alternatively, increasing pull rate extends the region of void formation to a larger radial core of the crystal and pushes the region of self-interstitials toward the periphery. For V ¼ 0:9 mm/ min, the self-interstitial rich region has been moved so close to the surface that diffusion to and recombination with at the surface of point defects has severely decreased the maximum in the intermediate CI concentration and self-interstitial clusters are not observed.
Fig. 6. Vacancy fields predicted at the pull rate values of 0.5 mm/min, 0.7 mm/min, and 0.9 mm/min.
106
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 8. Variation of V=GðRD¼0 Þ and V=GðRCV;max Þ with pull rate. The asymptotic result from Eq. (7) is shown for comparison.
This appearance of the peak suggests another definition for the location of the middle of the OSF-ring, as the location of the maximum in the residual vacancy profile after void nucleation; we define this value as RCV;max . This leads to ratio V ð8Þ ¼ K0max cm2 min1 K 1 GðRCV;max Þ
Fig. 7. Radial profiles of CV and CI computed (a) at z ¼ 2:8 cm (intermediate concentrations) and (b) near the top of the crystal (residual vacancy concentration) for several pull rates. The location of the OSF-ring corresponds to the region with high residual vacancy concentration in (b).
5. Residual vacancy concentrations, the osf-ring and almost perfect silicon An important feature of Fig. 7b is the appearance of the peak in the residual vacancy concentration at a particular value of the radius, which increases with pull rate.
as the more precise interpretation of the experiments of Dornberger et al. [14]. The expression Eq. (8) is plotted from our simulations on Fig. 8 and gives 1:50 103 4K0max 41:55 103 , compared to the value of K0 ¼ 1:43 103 computed from the asymptotic expression Eq. (7) and the range 1:28 103 4K0 41:35 103 extracted by plotting the value of V=GðRD¼0 Þ computed from the intermediate concentrations. The variation of K0 with both V=GðRD¼0 Þ and C=GðRCV;max Þ predicted from the simulations is partially a result of the realistic modeling of two-dimensional effects at the crystal surface introduced in the simulations. As shown in Fig. 6 and, certainly in the radial concentration profiles in Fig. 7a, the self-interstitial concentrations are affected by the presence of the crystal surface for pull rates as low as 0.7 mm/ min. Only in crystal growth of larger diameter crystals (800 diameter and above) would one expect
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
an appreciable range of V/G where K0 may be nearly constant. The values of K0 and K max are well within the 0 range of error inherent in estimating the interfacial temperature gradients G in the experiments. Moreover, the fact that the maximum in the residual vacancy concentration moves approximately with pull rate according to V=G suggests that this location can be recovered from an asymptotic analysis as well. However, such an analysis must include both the formation of voids and the effect of variations in the intermediate vacancy concentration. Analysis of the two-dimensional simulations also shows the fine-scale structure of the nearly defect-free region just outside of the OSF-ring. This structure is best seen by the plots in Fig. 9 of the radial structure of the total defects concentrations NV ðr; zÞ and NI ðr; zÞ computed near the top of the crystal for several pull rates. Also shown on these figures are the locations r ¼ RCV;max ðVÞ and r ¼ RD¼0 ðVÞ taken from the simulation for
107
V ¼ 0:6 mm/min. The region of low defect concentration is clearly visible near r ¼ Rfree ðVÞ and RCV;max 4Rfree 4RD¼0 . As with other markers for the radial structure, Rfree ðVÞvaries with V and according to V=G ðRfree Þ ¼ K0free cm2 =minK, where K0free is approximately equal to K0 . Again this result suggests that an asymptotic analysis can predict this dependence. The time-dependent microdefect dynamics simulations are useful tools to predict the spatial and temporal evolution of the microdefect structure with changes in operating conditions, such as the crystal pull rate V. The predicted evolution of the time-dependent dynamics of point defects and their clusters are shown in Fig. 10 for changing pull rate, but with a constant temperature field. The OSF-ring and IV boundaries separate as they approach the crystal surface with increasing pull rate. The number density of interstitial clusters decreases well as the IV neutral boundary comes close to the crystal surface and self-interstitials leave the crystal by diffusion and surface recombination. A lifetime map also is shown for (an experimental measure of the residual vacancy concentration in silicon) a horizontally sectioned CZ crystal in which a qualitatively similar growth rate schedule was used. The qualitative agreement in the two-dimensional structure of the boundary separating the void and self-interstitial regions is very good. More detailed comparison will have to await precise experiments and matching simulations.
6. Summary
Fig. 9. Radial profiles of total defect concentrations ðNV ðr; zÞ; NI ðr; zÞÞ computed at an axial position near the top of the crystal for pull rates (V) of 0.6 mm/min, 0.7 mm/min and 0.9 mm/min. Most significant is the prediction of an annular region of almost microdefect defect-free (of size greater than 50 nm) silicon as shown by the point (P) for V ¼ 0:7 mm/min, the locations of ðRD¼0 ; RCV;max ; Rfree Þ also are indicated for this pull rate.
The simulations reported here and in other references demonstrate the impact of quantitative analysis of defect dynamics on understanding microdefect formation during silicon crystal growth and subsequent processing. An enormous number of the patterns of microdefects seen during silicon crystal growth are semi-quantitatively predicted using self-consistent models of the equilibrium, transport, reaction and clustering of point defects. Quantitative agreement seems achievable, but is a strong function of accuracy of the thermophysical properties for point defects
108
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
Fig. 10. Simulation results for step change in crystal growth rate showing (a) vacancy density, and (b) lifetime map for crystal grown with similar step change in pull rate. Life time map of silicon crystal grown under similar conditions.
and of the level of detail incorporated into the simulations. The simulations reported here only described the interactions of point defects. Incorporation to include impurities, such as boron, is feasible [16]. More importantly, the extension of the analysis to include oxide precipitate formation has been accomplished [17]. When one-dimensional simulations including oxide precipitation are performed for temperature fields corresponding to the radial location of the maximum in residual vacancy concentration, r ¼ RCv;max , a measurable increase in oxygen precipitation is predicted during crystal cooling; these oxide precipitates are postulated to be the precursors of the defects in the OSF-ring. More interestingly, these defects are easily dissolved during a high temperature wafer annealing, leaving an almost defect-free (low values of NV and NI ) region in the crystal. Capitalizing on the appearance of the defect-free crystal region is difficult. To grow crystals with this
property across a wide portion of the radius would require designing a CZ hot-zone with an approximately radially uniform crystal temperature tailored so that defect clustering is avoided, independent of changes in pull rate and melt volume. This is a difficult challenge and it points to the need to couple the microdefect simulations to accurate macroscale global simulations of heat transfer and melt convection in commercial-scale CZ systems. The development of such simulations is underway [32]. The coupling through the temperature field in the crystal of these simulations with defect dynamics simulations will make the vision of simulation-aided defect engineering of large-scale CZ systems a reality [33].
Acknowledgements We thank MEMC Corporation, Shin-Etsu Hanotai Company and Wacker Siltronic for
R.A. Brown et al. / Journal of Crystal Growth 225 (2001) 97–109
financial support of this research. We are indebted to Erich Dornberger and Christian Hoess for the lifetime map of the crystal shown in Fig. 10.
References [1] A.J.R. de Kock, J. Wijgert, J. Crystal Growth 49 (1980) 718. [2] H. Yamagishi, I. Fusegawa, N. Fujimaki, M. Katayama, Semicond. Sci. Technol. 7 (1992) A135. [3] J. Ryuta, E. Morita, T. Tanaka, Y. Shimanuki, Jpn. J. Appl. Phys. 29 (1990) 1947. [4] M. Hourai, T. Nagashima, E. Kajita, S. Miki, T. Shigematsu, M. Okiu, J. Electrochem. Soc. 142 (1995) 3193. [5] T.S. Plaskett, Trans. Met. Soc. AIME, 233 (1965) 809. [6] H. Fo¨ll, U. Go¨sele, B.O. Kolbesen, J. Crystal Growth 52 (1981) 907. [7] T. Abe, H. Harada, MRS Symp. Proc. 14 (1983) 1. [8] T. Sinno, R.A. Brown, W. von Ammon, E. Dornberger, J. Electrochem. Soc. 145 (1998) 302. [9] V.V. Voronkov, J. Crystal Growth 59 (1982) 625. [10] R.A. Brown, D. Maroudas, T. Sinno, J. Crystal Growth 137 (1994) 12. [11] T. Sinno, R.A. Brown, J. Electrochem. Soc. 146 (1999) 2300. [12] T. Mori, Doctoral dissertation, Massachusetts Institute of Technology, Cambridge, MA, 2000. [13] T. Sinno, E. Dornberger, W. von Ammon, R.A. Brown, F. Dupret, Mater. Sci. Eng. 243 (2000) 1. [14] E. Dornberger, D. Gra¨f, M. Suhren, U. Lambert, P. Wagner, F. Dupret, W. von Ammon, J. Crystal Growth 180 (1997) 343. [15] T. Ueki, M. Itsumi, T. Takeda, Jpn. J. Appl. Phys. 38 (1999) 5695. [16] T. Sinno, H. Susanto, R.A. Brown, Appl. Phys. Lett. 75 (1999) 1544.
109
[17] Z. Wang, R.A. Brown, J. Crys. Growth (2000) submitted. [18] A.K. Kapila, Asymptotic Treatment of Chemically Reacting Systems, Pitman, London, 1983. [19] T. Mori, R.A. Brown, Chem. Eng. Sci. (2001) submitted. [20] T. Mori, Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, 2000. [21] J.S. Chang, G. Cooper, J. Comp. Phys. 6 (1970) 1. [22] B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods, Springer, New York, 1999. [23] K. Sueoka, M. Akatsuk, K. Nishihara, T. Yamamoto, S. Kobayashi, Mater. Sci. Forum 196–201 (1995) 1737. [24] M. Hasebe, Y. Takeoka, A. Shinoyama, S. Naito, J. Appl. Phys. 28 (1989) 1999. [25] V.V. Voronkov, R. Falster, J. Crystal Growth 194 (1998) 76. [26] V.V. Voronkov, R. Falster, J. Crystal Growth 204 (1999) 462. [27] M. Hourai, M. Nagashima, T. Kajita, E. Miki, S. Shigematsu, T. Okui, J. Electrochem Soc. 142 (1995) 3193. [28] K. Takano, K. Kitagawa, E. Iino, M. Kimura, H. Yamagishi, Mater. Sci. Forum 196 (1995) 1707. [29] R. Falster, V.V. Voronkov, J.C. Holzer, S. Markgraf, S. Macquaid, L. Mule’Stagno, In: H.R. Huff, H. Tsuya, U. Go¨sele, (Eds.), Semiconductor Silicon 1998, Electrochemical Society, Pennington, 1998. [30] M. Okiu, T. Tananka, T. Kanada, T. Ono, Ouyou Butsiri 66 (1997) 707. [31] T Saishouji, L. Nakamura, H. Nakajima, T. Yokoyama, F. Ishikawa, J. Tomioka, In: C.L. Claeys, P. RaiChoudhury, M. Watanabe, P. Stallhofer, H.J. Dawson, (Eds.), High Purity Silicon V, Electrochemical Society, Pennington, 1998. [32] A. Lipchin, R.A. Brown, J. Crystal Growth 216 (2000) 192–203. [33] R.A. Brown, Z. Wang, Y.C. Won, and A. Lipchin, J. Crystal Growth (2001) Invited paper, submitted.