- Email: [email protected]
Numerical modelling of the microscopic inhomogeneities during FZ silicon growth
Journal of Crystal Growth 198/199 (1999) 107—113
Numerical modelling of the microscopic inhomogeneities during FZ silicon growth A. Mu¨hlbauer *, A. Muiznieks , G. Raming , H. Riemann, A. Lu¨dge Institute for Electroheat, University of Hannover, Wilhelm-Busch-Strasse 4, D-30167 Hannover, Germany Institute of Crystal Growth, Rudower Chaussee 6, D-12489 Berlin, Germany
Abstract Transient axisymmetric numerical calculations of the hydrodynamic, temperature and solute concentration fields have been performed by means of FEM for the needle-eye FZ Silicon single-crystal growth process (diameter 4 ) to analyse the microscopic inhomogeneities. The rotation of the single crystal and feed rod, the buoyancy, Marangoni and electromagnetic (EM) forces in the melt are taken into account. Axisymmetric velocity oscillations caused by hydrodynamic instabilities are considered and calculated numerically. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: (1) the direct influence of the transient velocity field on the concentration field due to convective solute transport and (2) the influence of the oscillating temperature field on the local growth rate and as a consequence on the oscillating dopant segregation process at the growth interface. It is shown that for the considered experimental set-up the first mechanism dominates for the microscopic inhomogeneities. The calculated oscillations of the dopant concentration in the grown crystal (striations) are compared to spreading resistance measurements. 1999 Elsevier Science B.V. All rights reserved. Keywords: Microscopic inhomogeneities; Transient velocity; Numerical modelling; FZ silicon
1. Introduction The rapid progress in power semiconductor industry demands an increasing quality of the Silicon single crystals produced by FZ technique. The main goals are to ensure a good homogeneity of the macroscopic resistivity distributions and to avoid
* Corresponding author. Tel.: #49 511 762 2872; fax: #49 511 762 3275; e-mail: [email protected].
microscopic oscillations of resistivity (striations). Because of the complicated relation between the process parameters and the structure of striations in the grown crystal, mathematical modelling of the transient dopant incorporation is necessary. The mathematical model and the calculation method for the shape of the molten zone and hydrodynamics during FZ crystal growth is described in Ref. [1]. In Ref. [2] the calculation of the averaged-in-time dopant concentration fields and macroscopic resistivity distributions is given.
0022-0248/99/$ — see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 9 8 ) 0 1 1 5 8 - 0
108
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
Lu¨dge et al. [3] compare the results for macroscopic resistivity distributions calculated by the actual model to spreading resistance and 4-Point measurements. The numerical calculation of the transient dopant concentration fields in crystal growth is a new research field and is considered in a few publications. Ma and Walker [4] describe an asymptotic model for the unsteady transport of dopants during the growth of a semiconductor crystal from the melt with an externally applied magnetic field. The mass-diffusion boundary layers are considered and the dopant distribution in the entire crystal is predicted. Xiao [5] analysed the steady-state solute-concentration-driven flow, heat transfer and solute segregation during Czochralski growth of semiconductor compounds numerically using a finite-volume methodology. The effects of crystal rotation are also studied numerically. Jung and Mu¨ller [6] describe the numerical calculation of the transient, axisymmetric heat and solute transport including convection for a simplified vertical Bridgman process. Periodical oscillations were superimposed on the transient heater temperature profile. The amplitudes of the resulting oscillations of the growth rate and the dopant concentration (striations) in the growing crystals are compared with the predictions of analytical models.
In the present work the transient axisymmetric numerical calculations of the hydrodynamic, temperature and solute concentration fields were performed by means of FEM for the needle-eye FZ Silicon single crystal growth process (diameter 4 ) to analyse the microscopic inhomogeneities. In this content transient means that only short time scale transients relevant for the development of striations are considered for temperature, velocity and concentration in the melt phase. The large time-scale time-dependent effects, for example the change of shape of the solid and liquid phases during batch process are not taken into account. Since the time period considered in this paper (some minutes) is much smaller than the time constants due to batch process (about 5 h), both phenomenena can be decoupled. Fig. 1 shows an overview of the numerical calculations for the FZ process. The rotation of the crystal and feed rod, the buoyancy, Marangoni and EM forces in the melt are taken into account. Axisymmetric velocity oscillations caused by hydrodynamic instabilities of the analysed flow are considered and calculated numerically. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: (1) the direct influence of the transient velocity field on the concentration
Fig. 1. Finite-element mesh, calculated EM-, temperature, flow velocity fields and calculated radial resistivity distributions at different moments.
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
field due to convective solute transport and (2) the influence of the oscillating temperature field on the local growth rate and as a consequence on the oscillating dopant segregation process at the growth interface. It is shown, that for the considered experimental set-up the first mechanism dominates for the microscopic inhomogeneities. The calculated oscillations of the dopant concentration in the grown crystal (striations) are compared with the spreading resistance measurements.
2. Mathematical model The axisymmetric mathematical model for the calculation of the shape of the molten zone and the corresponding hydrodynamic (r, z, t) and temperature ¹(r, z, t) fields is described in Ref. [1]. The mathematical model and numerical procedure for the steady dopant concentration field are given in Ref. [2]. In case of transient dopant concentration field C(r, z, t) the following axisymmetric equation is used: *C #(v(r, z, t) ) )C"D*C, *t
(1)
where v(r, z, t) is the transient velocity field and D the dopant diffusivity. Because of possible crystal growth velocity oscillations v (r, t) the fol lowing boundary condition at the growth interface is used: *C D "v (r, t) ) (1!k ) ) C ) cos(a), *n
(2)
where a denotes the local angle between the growth interface and radial direction and k is the equilib rium segregation coefficient. n refers to the inner normal direction of the molten zone. In case of oscillating normal temperature derivative (*¹/*n)(r, t) in the melt at the growth interface, the oscillating growth velocity used in the Eq. (2) can be calculated by q (r)!j(*¹/*n)(r, t) v (r, t)" , oq cos(a)
(3)
109
where j denotes the melt heat conductivity, q the latent heat of fusion and o the density. The heat flux in the crystal at the growth interface q (r) is assumed to be time independent: *¹R q (r)"ov q cos (a)#j (r), *n
(4)
with (*¹R/*n)(r) as (*¹/*n)(r, t) averaged in time and v the pull velocity. Two mechanisms of the oscillat ing dopant incorporation in the crystal can be analysed by expressions (1) and (2). On the one hand there is a direct influence of the transient velocity field on the concentration field due to convective solute transport in Eq. (1). On the other hand, the oscillating temperature field influences the local growth rate (Eq. (3)) and as a consequence the oscillating dopant segregation process at the growth interface (boundary condition (2)). Possible origins of velocity oscillations in the melt during the FZ-process are the axisymmetric hydrodynamic instabilities of the rotating flow, the influence of the non-symmetric electromagnetic field of the HF-inductor and other reasons. In the present paper only the axisymmetric hydrodynamic instabilities for the experimental set-up depicted in Ref. [1] are investigated numerically and the corresponding transient dopant concentration fields are calculated.
3. Numerical method The numerical method for the calculation of the shape of the molten zone and the corresponding hydrodynamic and temperature fields is described in Ref. [1]. In Ref. [2] the numerical method and the used finite-element mesh for the concentration field with distinct boundary layer is presented. Because of a thin solute boundary layer in the vicinity of the growth interface an essentially finer finite-element grid than in case of hydrodynamic calculation is used in order to deal with the highconcentration gradients more accurately. The calculated transient velocity field is interpolated to obtain the velocity values on this finer mesh near the growth interface. The transient calculations are performed for periods of about 100 s and the fields
110
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
at some thousand different moments are stored during hydrodynamic and dopant concentration field calculation. The time-dependent dopant concentration fields are averaged numerically to yield the macroscopic resistivity distribution in the single crystal. For the microscopic resistivity oscillations the concentration distribution at the growth interface at different moments and corresponding axial positions is used.
4. Results of calculations and comparison with experiment The presented calculations are carried out for the FZ Silicon growth laboratory system at Institute of Crystal Growth, Berlin. The used physical properties of Si and main growth parameters are described in Table 1. For more details see Refs. [1,2]. First, the shape of the molten zone is calculated and calculations of the transient coupled hydrodynamic v(r, z, t) and thermal ¹(r, z, t) fields are carried out. Starting with zero fields, after a time period of approximately 100 s a fully developed oscillating hydrodynamic and thermal solution is reached. For the following 100 s the oscillating velocity v(r, z, t) and temperature fields ¹(r, z, t) are stored and used for concentration calculations. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: First, the direct influence of the transient velocity field v(r, z, t) on the concentration field due to convective solute transport is investigated. Therefore, it is assumed that the local growth rate is time independent and equals the pull velocity: v (r, t)"v . Fig. 2 shows the calculated hydrodynamic fields (a) and dopant concentration fields (b) at five subsequent moments for the fully developed hydrodynamic situation in case of a constant local growth rate and crystal rotation rate 5 rpm. During the crystallization process, the increased dopant concentration (k (1) in the melt is obtained directly at the growth interface in the diffusion layer. This increased dopant concentration is affected by the unsteady bulk flow in the melt with a lower dopant concentration. The oscillating concentration at the growth interface is used
Table 1 Physical properties and process parameters Diffusivity Melt heat conductivity Segregation coefficient Solid conductivity Liquid conductivity Latent heat of fusion Expansion coefficient Kinematic viscosity Melting temperature Heat capacity Density Growth velocity HF-inductor current Current frequency Crystal rotation Feed rotation Crystal radius
3.4;10\ m/s 67 W/mK 0.35 5;10 1/)m 1.2;10 1/)m 1.8;10 J/kg 1.4;10\ 1/K 3.4;10\ m/s 1685 K 1000 J/(kgK) 2530 kg/m 3.4 mm/min 990 A 2.8 MHz 5, 10 and 15 rpm !15 rpm 50 mm
to calculate the normalized resistivity in the single crystal as the reciprocal value of the concentration. By conversion of time to axial coordinate in the crystal the resistivity oscillations in axial direction (striations) are calculated (Fig. 3a). The second mechanism leading to resistivity striations is determined by the influence of the oscillating temperature field ¹(r, z, t) on the local growth rate v (r, t) and as a consequence on the oscillating dopant segregation process at the growth interface. Fig. 4a shows the variation in time of the local growth rate v (r, t) in case of crystal rotation rate 5 rpm calculated from ¹(r, z, t) by expression (3). The corresponding calculation of C(r, z, t) is carried out by using a time-independent velocity field (v(r, z, t) averaged over 50 s) for convective solute transport. The resistivity oscillations in the grown crystal calculated from oscillating concentration at the growth interface are shown in Fig. 4b. It can be seen that for the considered set-up the resistivity oscillations caused by the time-dependent local growth velocity are much smaller than that caused directly by velocity oscillations (Fig. 3a). Therefore, the influence of the time-dependent local growth velocity has been neglected for further considerations. Fig. 3 also shows the comparison between spreading resistance measurements and calculated
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
111
Fig. 2. Flow velocity fields (a, streamlines with increment 5;10\ m/s) and normalized dopant concentration fields (b, isolines with increment 0.05, molten material from feed rod has concentration 1) at subsequent moments.
112
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
Fig. 3. Comparison of calculated (a) and measured (b) microscopic resistivity profiles in axial direction near the crystal center, at half radius and near the rim.
microscopic resistivity distributions in axial direction. For comparison three representative positions near the crystal center, at half radius and near the rim have been chosen for both experiment and calculation. The axial distance of 6 mm in grown crystal corresponds to a calculation time of about 106 s at a crystal growth rate of 3.4 mm/min. In the central region of the crystal the amplitude of measured resistivity oscillations (Fig. 3b) is much higher than the oscillation amplitude at the periphery. The absolute level of resistivity in the central region is about two times higher than at the periphery. This is confirmed by the numerical model (Fig. 3).
Fig. 4. Variation in time of the local growth rate (a) and corresponding resistivity oscillations (b) in the single crystal in axial direction.
The frequency of resistivity oscillations for measurement and calculation is in agreement, too. However, as Fig. 3 shows, the amplitudes of calculated oscillations are smaller than the measured ones. Especially in peripheral regions the measurements show an influence of the non-symmetric electromagnetic field due to the current suppliers of the inductor that cause the so-called rotational striations. In the central part, the calculated oscillations could be smaller than measured ones because the axisymmetric model does not cover the evolution of three-dimensional oscillations. To overcome the mentioned limitations of
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
113
lation show that at higher rotation rates the minimum of the resistivity distribution is shifted to the periphery of the crystal.
5. Conclusions The axisymmetric mathematical model presented here covers the whole FZ Silicon crystal growth process from calculation of the interface shape and melt flow up to the dopant concentration field in the melt in dependence on process parameters. The calculated transient concentration in the molten zone at the growth interface is used for the calculation of macroscopic and microscopic resistivity distributions in the single crystal. A comparison between the calculated macroscopic and microscopic resistivity distributions and measured ones led to a good agreement. As known from practice, the resistivity oscillations in the central region of the single crystal are more intense than in the outer region of the crystal. This has been confirmed by the numerical calculations. The calculated examples show that for hydrodynamic conditions of the considered FZ-equipment the concentration oscillations determined directly by the fluctuating velocity field overcome concentration oscillations caused by the oscillating growth velocity. Fig. 5. Comparison of measured and calculated macroscopic resistivity profiles at crystal rotation rates 10 rpm (a) and 15 rpm (b).
the actual 2D-model, 3D-calculations of the EM-, HD- and dopant concentration fields are in progress. The calculated unsteady concentration fields can be used to yield the averaged-in-time concentration distribution at the growth interface. The resulting normalized macroscopic resistivity distribution is calculated as the reciprocal value of the averaged dopant concentration. Fig. 5 shows the comparison of measured (4-Point measurements) and calculated macroscopic resistivity distributions at rotation rates 10 and 15 rpm. The shape of the profile, the position of the minimum and the levels coincide very well. Both, experiment and calcu-
Acknowledgements The authors are grateful for the support received from Wacker-Siltronic AG, Burghausen, Germany.
References [1] A. Mu¨hlbauer, A. Muiznieks, J. Virbulis, A. Lu¨dge, H. Riemann, J. Crystal Growth 151 (1995) 66. [2] A. Mu¨hlbauer, A. Muiznieks, J. Virbulis, J. Crystal Growth 180 (1997) 372. [3] A. Lu¨dge, H. Riemann, W. Schro¨der, A. Mu¨hlbauer, A. Muiznieks, G. Raming, 8th Int. Symp. on Silicon Materials Science and Technology, 3—8 May 1998, San Diego. [4] N. Ma, J.S. Walker, J. Crystal Growth 172 (1997) 124. [5] Q. Xiao, J. Crystal Growth 174 (1997) 7. [6] T. Jung, G. Mu¨ller, J. Crystal Growth 171 (1997) 373.
Numerical modelling of the microscopic inhomogeneities during FZ silicon growth A. Mu¨hlbauer *, A. Muiznieks , G. Raming , H. Riemann, A. Lu¨dge Institute for Electroheat, University of Hannover, Wilhelm-Busch-Strasse 4, D-30167 Hannover, Germany Institute of Crystal Growth, Rudower Chaussee 6, D-12489 Berlin, Germany
Abstract Transient axisymmetric numerical calculations of the hydrodynamic, temperature and solute concentration fields have been performed by means of FEM for the needle-eye FZ Silicon single-crystal growth process (diameter 4 ) to analyse the microscopic inhomogeneities. The rotation of the single crystal and feed rod, the buoyancy, Marangoni and electromagnetic (EM) forces in the melt are taken into account. Axisymmetric velocity oscillations caused by hydrodynamic instabilities are considered and calculated numerically. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: (1) the direct influence of the transient velocity field on the concentration field due to convective solute transport and (2) the influence of the oscillating temperature field on the local growth rate and as a consequence on the oscillating dopant segregation process at the growth interface. It is shown that for the considered experimental set-up the first mechanism dominates for the microscopic inhomogeneities. The calculated oscillations of the dopant concentration in the grown crystal (striations) are compared to spreading resistance measurements. 1999 Elsevier Science B.V. All rights reserved. Keywords: Microscopic inhomogeneities; Transient velocity; Numerical modelling; FZ silicon
1. Introduction The rapid progress in power semiconductor industry demands an increasing quality of the Silicon single crystals produced by FZ technique. The main goals are to ensure a good homogeneity of the macroscopic resistivity distributions and to avoid
* Corresponding author. Tel.: #49 511 762 2872; fax: #49 511 762 3275; e-mail: [email protected].
microscopic oscillations of resistivity (striations). Because of the complicated relation between the process parameters and the structure of striations in the grown crystal, mathematical modelling of the transient dopant incorporation is necessary. The mathematical model and the calculation method for the shape of the molten zone and hydrodynamics during FZ crystal growth is described in Ref. [1]. In Ref. [2] the calculation of the averaged-in-time dopant concentration fields and macroscopic resistivity distributions is given.
0022-0248/99/$ — see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 9 8 ) 0 1 1 5 8 - 0
108
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
Lu¨dge et al. [3] compare the results for macroscopic resistivity distributions calculated by the actual model to spreading resistance and 4-Point measurements. The numerical calculation of the transient dopant concentration fields in crystal growth is a new research field and is considered in a few publications. Ma and Walker [4] describe an asymptotic model for the unsteady transport of dopants during the growth of a semiconductor crystal from the melt with an externally applied magnetic field. The mass-diffusion boundary layers are considered and the dopant distribution in the entire crystal is predicted. Xiao [5] analysed the steady-state solute-concentration-driven flow, heat transfer and solute segregation during Czochralski growth of semiconductor compounds numerically using a finite-volume methodology. The effects of crystal rotation are also studied numerically. Jung and Mu¨ller [6] describe the numerical calculation of the transient, axisymmetric heat and solute transport including convection for a simplified vertical Bridgman process. Periodical oscillations were superimposed on the transient heater temperature profile. The amplitudes of the resulting oscillations of the growth rate and the dopant concentration (striations) in the growing crystals are compared with the predictions of analytical models.
In the present work the transient axisymmetric numerical calculations of the hydrodynamic, temperature and solute concentration fields were performed by means of FEM for the needle-eye FZ Silicon single crystal growth process (diameter 4 ) to analyse the microscopic inhomogeneities. In this content transient means that only short time scale transients relevant for the development of striations are considered for temperature, velocity and concentration in the melt phase. The large time-scale time-dependent effects, for example the change of shape of the solid and liquid phases during batch process are not taken into account. Since the time period considered in this paper (some minutes) is much smaller than the time constants due to batch process (about 5 h), both phenomenena can be decoupled. Fig. 1 shows an overview of the numerical calculations for the FZ process. The rotation of the crystal and feed rod, the buoyancy, Marangoni and EM forces in the melt are taken into account. Axisymmetric velocity oscillations caused by hydrodynamic instabilities of the analysed flow are considered and calculated numerically. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: (1) the direct influence of the transient velocity field on the concentration
Fig. 1. Finite-element mesh, calculated EM-, temperature, flow velocity fields and calculated radial resistivity distributions at different moments.
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
field due to convective solute transport and (2) the influence of the oscillating temperature field on the local growth rate and as a consequence on the oscillating dopant segregation process at the growth interface. It is shown, that for the considered experimental set-up the first mechanism dominates for the microscopic inhomogeneities. The calculated oscillations of the dopant concentration in the grown crystal (striations) are compared with the spreading resistance measurements.
2. Mathematical model The axisymmetric mathematical model for the calculation of the shape of the molten zone and the corresponding hydrodynamic (r, z, t) and temperature ¹(r, z, t) fields is described in Ref. [1]. The mathematical model and numerical procedure for the steady dopant concentration field are given in Ref. [2]. In case of transient dopant concentration field C(r, z, t) the following axisymmetric equation is used: *C #(v(r, z, t) ) )C"D*C, *t
(1)
where v(r, z, t) is the transient velocity field and D the dopant diffusivity. Because of possible crystal growth velocity oscillations v (r, t) the fol lowing boundary condition at the growth interface is used: *C D "v (r, t) ) (1!k ) ) C ) cos(a), *n
(2)
where a denotes the local angle between the growth interface and radial direction and k is the equilib rium segregation coefficient. n refers to the inner normal direction of the molten zone. In case of oscillating normal temperature derivative (*¹/*n)(r, t) in the melt at the growth interface, the oscillating growth velocity used in the Eq. (2) can be calculated by q (r)!j(*¹/*n)(r, t) v (r, t)" , oq cos(a)
(3)
109
where j denotes the melt heat conductivity, q the latent heat of fusion and o the density. The heat flux in the crystal at the growth interface q (r) is assumed to be time independent: *¹R q (r)"ov q cos (a)#j (r), *n
(4)
with (*¹R/*n)(r) as (*¹/*n)(r, t) averaged in time and v the pull velocity. Two mechanisms of the oscillat ing dopant incorporation in the crystal can be analysed by expressions (1) and (2). On the one hand there is a direct influence of the transient velocity field on the concentration field due to convective solute transport in Eq. (1). On the other hand, the oscillating temperature field influences the local growth rate (Eq. (3)) and as a consequence the oscillating dopant segregation process at the growth interface (boundary condition (2)). Possible origins of velocity oscillations in the melt during the FZ-process are the axisymmetric hydrodynamic instabilities of the rotating flow, the influence of the non-symmetric electromagnetic field of the HF-inductor and other reasons. In the present paper only the axisymmetric hydrodynamic instabilities for the experimental set-up depicted in Ref. [1] are investigated numerically and the corresponding transient dopant concentration fields are calculated.
3. Numerical method The numerical method for the calculation of the shape of the molten zone and the corresponding hydrodynamic and temperature fields is described in Ref. [1]. In Ref. [2] the numerical method and the used finite-element mesh for the concentration field with distinct boundary layer is presented. Because of a thin solute boundary layer in the vicinity of the growth interface an essentially finer finite-element grid than in case of hydrodynamic calculation is used in order to deal with the highconcentration gradients more accurately. The calculated transient velocity field is interpolated to obtain the velocity values on this finer mesh near the growth interface. The transient calculations are performed for periods of about 100 s and the fields
110
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
at some thousand different moments are stored during hydrodynamic and dopant concentration field calculation. The time-dependent dopant concentration fields are averaged numerically to yield the macroscopic resistivity distribution in the single crystal. For the microscopic resistivity oscillations the concentration distribution at the growth interface at different moments and corresponding axial positions is used.
4. Results of calculations and comparison with experiment The presented calculations are carried out for the FZ Silicon growth laboratory system at Institute of Crystal Growth, Berlin. The used physical properties of Si and main growth parameters are described in Table 1. For more details see Refs. [1,2]. First, the shape of the molten zone is calculated and calculations of the transient coupled hydrodynamic v(r, z, t) and thermal ¹(r, z, t) fields are carried out. Starting with zero fields, after a time period of approximately 100 s a fully developed oscillating hydrodynamic and thermal solution is reached. For the following 100 s the oscillating velocity v(r, z, t) and temperature fields ¹(r, z, t) are stored and used for concentration calculations. Two mechanisms of the oscillating dopant incorporation in the crystal are investigated: First, the direct influence of the transient velocity field v(r, z, t) on the concentration field due to convective solute transport is investigated. Therefore, it is assumed that the local growth rate is time independent and equals the pull velocity: v (r, t)"v . Fig. 2 shows the calculated hydrodynamic fields (a) and dopant concentration fields (b) at five subsequent moments for the fully developed hydrodynamic situation in case of a constant local growth rate and crystal rotation rate 5 rpm. During the crystallization process, the increased dopant concentration (k (1) in the melt is obtained directly at the growth interface in the diffusion layer. This increased dopant concentration is affected by the unsteady bulk flow in the melt with a lower dopant concentration. The oscillating concentration at the growth interface is used
Table 1 Physical properties and process parameters Diffusivity Melt heat conductivity Segregation coefficient Solid conductivity Liquid conductivity Latent heat of fusion Expansion coefficient Kinematic viscosity Melting temperature Heat capacity Density Growth velocity HF-inductor current Current frequency Crystal rotation Feed rotation Crystal radius
3.4;10\ m/s 67 W/mK 0.35 5;10 1/)m 1.2;10 1/)m 1.8;10 J/kg 1.4;10\ 1/K 3.4;10\ m/s 1685 K 1000 J/(kgK) 2530 kg/m 3.4 mm/min 990 A 2.8 MHz 5, 10 and 15 rpm !15 rpm 50 mm
to calculate the normalized resistivity in the single crystal as the reciprocal value of the concentration. By conversion of time to axial coordinate in the crystal the resistivity oscillations in axial direction (striations) are calculated (Fig. 3a). The second mechanism leading to resistivity striations is determined by the influence of the oscillating temperature field ¹(r, z, t) on the local growth rate v (r, t) and as a consequence on the oscillating dopant segregation process at the growth interface. Fig. 4a shows the variation in time of the local growth rate v (r, t) in case of crystal rotation rate 5 rpm calculated from ¹(r, z, t) by expression (3). The corresponding calculation of C(r, z, t) is carried out by using a time-independent velocity field (v(r, z, t) averaged over 50 s) for convective solute transport. The resistivity oscillations in the grown crystal calculated from oscillating concentration at the growth interface are shown in Fig. 4b. It can be seen that for the considered set-up the resistivity oscillations caused by the time-dependent local growth velocity are much smaller than that caused directly by velocity oscillations (Fig. 3a). Therefore, the influence of the time-dependent local growth velocity has been neglected for further considerations. Fig. 3 also shows the comparison between spreading resistance measurements and calculated
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
111
Fig. 2. Flow velocity fields (a, streamlines with increment 5;10\ m/s) and normalized dopant concentration fields (b, isolines with increment 0.05, molten material from feed rod has concentration 1) at subsequent moments.
112
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
Fig. 3. Comparison of calculated (a) and measured (b) microscopic resistivity profiles in axial direction near the crystal center, at half radius and near the rim.
microscopic resistivity distributions in axial direction. For comparison three representative positions near the crystal center, at half radius and near the rim have been chosen for both experiment and calculation. The axial distance of 6 mm in grown crystal corresponds to a calculation time of about 106 s at a crystal growth rate of 3.4 mm/min. In the central region of the crystal the amplitude of measured resistivity oscillations (Fig. 3b) is much higher than the oscillation amplitude at the periphery. The absolute level of resistivity in the central region is about two times higher than at the periphery. This is confirmed by the numerical model (Fig. 3).
Fig. 4. Variation in time of the local growth rate (a) and corresponding resistivity oscillations (b) in the single crystal in axial direction.
The frequency of resistivity oscillations for measurement and calculation is in agreement, too. However, as Fig. 3 shows, the amplitudes of calculated oscillations are smaller than the measured ones. Especially in peripheral regions the measurements show an influence of the non-symmetric electromagnetic field due to the current suppliers of the inductor that cause the so-called rotational striations. In the central part, the calculated oscillations could be smaller than measured ones because the axisymmetric model does not cover the evolution of three-dimensional oscillations. To overcome the mentioned limitations of
A. Mu¨ hlbauer et al. / Journal of Crystal Growth 198/199 (1999) 107–113
113
lation show that at higher rotation rates the minimum of the resistivity distribution is shifted to the periphery of the crystal.
5. Conclusions The axisymmetric mathematical model presented here covers the whole FZ Silicon crystal growth process from calculation of the interface shape and melt flow up to the dopant concentration field in the melt in dependence on process parameters. The calculated transient concentration in the molten zone at the growth interface is used for the calculation of macroscopic and microscopic resistivity distributions in the single crystal. A comparison between the calculated macroscopic and microscopic resistivity distributions and measured ones led to a good agreement. As known from practice, the resistivity oscillations in the central region of the single crystal are more intense than in the outer region of the crystal. This has been confirmed by the numerical calculations. The calculated examples show that for hydrodynamic conditions of the considered FZ-equipment the concentration oscillations determined directly by the fluctuating velocity field overcome concentration oscillations caused by the oscillating growth velocity. Fig. 5. Comparison of measured and calculated macroscopic resistivity profiles at crystal rotation rates 10 rpm (a) and 15 rpm (b).
the actual 2D-model, 3D-calculations of the EM-, HD- and dopant concentration fields are in progress. The calculated unsteady concentration fields can be used to yield the averaged-in-time concentration distribution at the growth interface. The resulting normalized macroscopic resistivity distribution is calculated as the reciprocal value of the averaged dopant concentration. Fig. 5 shows the comparison of measured (4-Point measurements) and calculated macroscopic resistivity distributions at rotation rates 10 and 15 rpm. The shape of the profile, the position of the minimum and the levels coincide very well. Both, experiment and calcu-
Acknowledgements The authors are grateful for the support received from Wacker-Siltronic AG, Burghausen, Germany.
References [1] A. Mu¨hlbauer, A. Muiznieks, J. Virbulis, A. Lu¨dge, H. Riemann, J. Crystal Growth 151 (1995) 66. [2] A. Mu¨hlbauer, A. Muiznieks, J. Virbulis, J. Crystal Growth 180 (1997) 372. [3] A. Lu¨dge, H. Riemann, W. Schro¨der, A. Mu¨hlbauer, A. Muiznieks, G. Raming, 8th Int. Symp. on Silicon Materials Science and Technology, 3—8 May 1998, San Diego. [4] N. Ma, J.S. Walker, J. Crystal Growth 172 (1997) 124. [5] Q. Xiao, J. Crystal Growth 174 (1997) 7. [6] T. Jung, G. Mu¨ller, J. Crystal Growth 171 (1997) 373.