The thermal analysis of the Czochralski crystal growth process
is very difficult to accurately model, due to the interaction of heat conduction,
convection-diffusion, thermal radiation, fluid flow, and other transport
phenomena.
In collaboration with the Consortium for Crystal Growth funded by DARP/AFOSR
a highly innovative, general-purpose computer code for phase-change and
free-surface problems, utilizing multi-zone adaptive grid generation and
a curvilinear finite volume scheme, is linked to a volume and surface thermal
radiation code to evaluate the temperature distribution within a Czochralski
crystal furnace for several cases. The results of the present study
show that the temperature distribution within the crystal changes with
crystal's optical thickness as well as scattering albedo and interface
location.
BACKGROUND
Control of crystal growth processes is limited by the degree to which the complex non-linear transport phenomena involved is understood. The heat transfer environment, established by conductive, convective and radiative modes of energy transfer, must yield a well-controlled temperature field favorable to conditions for the growth of low-defect, high-quality crystals. In many optical crystals, such as yttrium aluminum garnet (YAG), gadolinium gallium garnet (GGG), and sapphire (Al2O3), radiative exchange within the crystalline medium, termed internal radiation, is considered important (Brandon and Derby, 1991). Cockayne et al. (1969) and Kvapil et al. (1981) report evidence correlating crystal/melt interface shape with internal radiative transport during YAG growth. Moreover, Yuferev and Vasil'ev (1987) have shown that radiative transfer within sapphire ribbon can lead to instability of the crystallization front. Due to the inherent complexities associated with radiative transfer calculations, such as the long distance nature of radiation (solid angle integration) and functional dependence of radiative properties, radiative analyses of systems with semi-transparent crystals are very often simplified. Consequently, simplification or neglect of detailed radiative phenomena may inaccurately predict the temperature field, interface shape/location, and thermal stress distribution within the growth furnace.
Semiconductor and optical single crystals are commonly grown by the Czochralski (CZ) method. The CZ process is performed in a crystal growing apparatus, illustrated in Fig. 1, which consists of a quartz crucible, graphite susceptor, and rotating graphite pedestal. In CZ processes, a seed crystal is dipped into a molten pool, whereafter the thermal environment is dynamically adjusted so that a crystal grows from the seed as it is pulled from the melt. A crystal grown during the CZ process exchanges thermal radiation (medium is non-participating) with the melt surface, crucible wall, and puller interior surfaces, all of which are at different temperatures and exhibit different surface emissivities. Various approaches have been utilized to model surface radiation phenomena in CZ enclosures, including isothermal treatment of surfaces excluding the crystal (Rea, 1981), establishment of radiation heat transfer coefficients to linearize boundary conditions (Williams and Reuser, 1983; Motakef and Witt, 1987), and use of generalized view factor programs, such as FACET (Shapiro, 1983; Atherton, 1987). In some instances, it is plausible to approximate the radiative heat flux through spectral, semi-transparent materials by dividing the wavelength spectrum into two bands, for which the radiative medium is either totally transparent (zero absorption coefficient) or completely opaque (very large absorption coefficient). This methodology was used by Crochet et al. (1989) to model radiative exchanges within the quartz crucible lining in vertical Bridgman crystal growth processes. In liquid-encapsulated Czochralski (LEC) crystal growth, a highly viscous oxide layer, usually boric oxide B2O3, is layered over the melt to suppress the loss of volatile components in high pressure systems such as indium phosphate (InP) and gallium arsenide (GaAs). The B2O3 encapsulant is spectral/semi-transparent and has been modeled using the two-band approximation in global simulations of LEC GaAs growth by Dupret et al. (1989, 1990) and Thomas et al. (1989). In fact, many types of crystalline materials demonstrate the indicated two-band optical behavior (Modest, 1993). In an effort to approximate the effects of internal radiative transport within yttrium aluminum garnet (YAG) and gadolinium gallium garnet (GGG), Xiao and Derby (1994) considered the interior of the crystal as a radiative enclosure, bound by vanishingly thin opaque surfaces. This simplified treatment of opaque surfaces is valid for real systems, where surfaces are often coated with impurities deposited by vapor-phase transport (1993).

Traditionally, models of crystal growth have simplified the effects of internal radiative transfer on the growth process, as in the two-band approximation. This approach is justified for the treatment of most semiconductor growth systems, such as Si and GaAs. However, the optical properties of semi-transparent crystals indicate that internal radiation within the crystal should not be neglected (Brandon and Derby, 1991). Analyses of internal radiative transport within crystal growth systems have in the past been limited to idealized, one-dimensional conduction and radiation studies (O'Hara et al., 1968; Abrams and Viskanta, 1974). On a more recent note, several investigators have rigorously modeled internal radiation within the CZ system. Matsushima and Viskanta (1990) used the Spherical Harmonics (P1) method to approximate radiative heat transfer in a semi-transparent crystal and the melt. Tsukada et al. (1995) used a similar procedure to study the effect of the absorption coefficient on the crystal/melt interface shape for lithium niobate (LiNbO3). The P-1 method reduces the integral equations of radiative intensity for a participating media to a set of simultaneous partial differential equations, thereby facilitating the solution of the constituent governing equations. However, this differential approximation is based on the assumption of an optically thick medium and can yield results serious in error for optically thin regions. Brandon and Derby (1991, 1992) employed a Galerkin finite element strategy and three/nine point Gaussian integration routine to directly integrate the radiative transport equation at each Gaussian point over the solid angle. The computations were performed on a Cray X-MP supercomputer and are extremely arduous/expensive, since three-dimensional information is needed to perform the solid angle integration. Overall, numerical simulations of crystal growth processes require comprehensive treatment of the thermophysical phenomena involved. Since crystal growth is a transient phase-change process with moving, irregularly-shaped boundaries, radiative transfer models should be computationally flexible so as to deal with the fluctuating shape of the crystal/melt interface, formation of meniscii at free-surfaces, and blockage of view by obstructing bodies. Furthermore, radiative calculations are performed at each iteration/time step, making computational efficiency a key issue.
A surface radiation module (Nunes et al., 1996) and volume radiation algorithm, based on the Discrete Exchange Factor (DEF) method , for radiative analyses of irregularly-shaped axisymmetric enclosures with participating media (Nunes and Naraghi, 1997) is combined with a multi-zone adaptive transport code, MASTRAPP2d (Zhang and Prasad, 1995b), to simulate CZ crystal growth. MASTRAPP2d utilizes a multi-zone adaptive grid-generation technique (MAGG) for discretizing the physical domain and a curvilinear finite volume (CFV) technique for solving the governing equations. The MAGG algorithm distributes grid nodes adaptively in response to the development of the solution and according to the evolution of the problem domain. Further details on the numerical methodology of MASTRAPP2d is provided by Zhang and Prasad (1995a). A schematic of the CZ system used in simulating the growth process is shown in Fig. 2. For the system at hand, two separate radiative enclosures are considered, each bound by thin opaque surfaces (Xiao and Derby, 1994). The first radiative enclosure is comprised of the outer surface of the crystal body, the free-surface of the melt, the exposed portion of the crucible wall, and an enclosing top surface. The components of the second radiative enclosure include the crystal cap, crystal body, and the crystal/melt interface. All surfaces are subdivided into a finite number of isothermal differential ring elements for which radiative heat fluxes are ultimately computed. The surface radiation algorithm is incorporated as a subroutine in MASTRAPP2d to evaluate the radiative heat fluxes at the surfaces of the crystal body, free-surface of the melt, exposed portion of the crucible wall, and enclosing top surface (radiative enclosure #1 in Fig. 2). Internal radiative effects within the crystal (radiative enclosure #2 in Fig. 2) are modeled using the volume radiation algorithm. The DEF method, based upon a point-to-point approach for radiative analyses of enclosures has proven to be a computationally efficient and accurate method (Naraghi et al., 1988). Direct exchange factors between differential ring element pairs are computed by generalizing Modest's (1988, 1993) model for view factors between differential ring elements on concentric bodies to give the appropriate DEF expressions. Total exchange factors are subsequently computed to allow for multiple reflections and scattering of radiation by surfaces/media. Blockage effects produced by inner and/or outer bodies are accounted for in the present formulation. The codes are used for purposes of theoretical investigation of the effects of extinction coefficient, scattering albedo, and boundary conditions upon the furnace temperature distribution and crystal/melt interface movement for the CZ growth of semi-transparent crystals. A complete description of the mathematical model is given by Nunes and Naraghi (1997 and 1998)

RESULTS AND DISCUSSION
The surface/volume radiation codes are linked to an axisymmetric adaptive transport scheme based on MAGG and CFV procedures to study the effect of thermal radiative transport on the CZ crystal growth of silicon. A parametric study of the effects of extinction coefficient, scattering albedo, and boundary condition type at the crystal cap is presented. All numerical simulations were performed on a Silicon Graphics workstation, for 80 time steps, requiring approximately 8 hours of computational time for a run. Several assumptions are made for these computations: (1) velocity at the free surface of the melt due to a change in melt height is negligible; (2) the thermophysical and radiative properties of the crystal and melt are invariant; (3) the crystal pull rate is varied to produce a crystal of constant diameter; (4) all surfaces are diffuse and opaque; (5) the melt region is opaque; (6) index of refraction of crystalline medium is that of silicon (3.42); (7) the emissivity of the crystal, free-surface of the melt, crucible wall, enclosing top surface, and crystal/melt interface are 0.7, 0.3, 0.8, 0.5 and 0.5, respectively (Sackinger et al., 1989); (8) internal radiative transport within the semitransparent crystal is spectral and exists for a band with limits of 0 - 4.75 micro-meters. All non-radiative transport parameters are fixed for the analyses, and are given in (Nunes and Naraghi, 1998). A sample mesh layout is provided in Fig. 3.

For the purposes of theoretical investigation, consideration is first given to a short crystal with prescribed crystal cap temperature. Figures 4-7 give the dimensionless temperature contours for values of extinction coefficient within the crystalline medium of Kt= 0.1, 1, 5, and 20m-1, respectively. The scattering albedo for these cases is set to zero. Due to the imposed temperatures at the crystal cap and crystal/melt interface (represented by the "0" contour), relatively large amounts of radiant energy are transferred form the interface and lower portion of the crystal to the cap for low values of extinction coefficient (Figs. (4) and (5)). Consequently, large radial temperature gradients are evident towards the foot of the crystal, especially by the surface. As radiant emission from the interface is allowed to decrease (by increasing the extinction coefficient), as shown in Figs. (6) and (7), the interface behaves counter-intuitively and shifts downwards. Interestingly, as less energy is transferred from the bottom of the crystal to the cap, the radial temperature gradients diminish. Consequently, radial conductive heat fluxes in the lower section of the crystal are decreased, which act to deflect the interface downwards. Additionally, as the extinction coefficient of the medium increases, external radiative effects (emission from free-surface, crucible wall, and top surface) start to dominate over internal radiative transport.




A more realistic top cap boundary condition is the conductive-radiative condition the same as that of the cylindrical crystal surface given. Note, the crystal for this case is longer than the previous cases. The results for Kt= 0 (transparent crystal), Kt= 1 and 10m-1, with zero scattering, are shown in Figs. 8-10. For low values of extinction coefficient (Figs. (8) and (9)), as in the previous case, radial temperature gradients develop at the lower section of the crystal, with evidence of absorption at the top of the crystal. As the extinction coefficient is increased further, less energy penetrates through the interior of the crystal from the interface and external radiative effects become noticeable. The crystal/melt interface, however, shows relatively little movement, due to the low axial temperature gradient (approximately 40oC between the interface and cap) allowed by the conduction-radiation condition at the cap.



Next, the effects of varying the scattering albedo on the temperature distribution within the crystal are examined. Figures (11) and (12) gives the dimensionless temperature distribution for Kt= 1 and 10m-1, with a scattering albedo of 0.5. There is very little difference between these contours and those for a purely absorbing/emitting medium. However, if the scattering albedo is increased to 1.0 and the medium is purely scattering (Figs. (13) and (14)), the temperature distribution within the crystal is markedly different. Interestingly, radiant emission at the foot of the crystal for the latter cases in notably higher, with very low radial gradients towards the middle and upper sections. Again, the overall axial temperature gradient for these cases are very low.




In crystals, the radiative properties are highly temperature dependent.
The temperature and spatial variations make solving the radiative transport
equations extremely difficult. Not only are the governing equations
highly nonlinear requiring an increased number of iterations, but also
calculations of the transmittance is complicated by the fact
that the extinction coefficient varies along the optical pathlength.
Preliminary work has been done by Neghabat and Naraghi (1998) on a rectangular
non-homogeneous media. Our future work will be to extend the model to irregular
and axisymmetric configurations and finally link the model to the MASTRAPP.
REFERENCES
Abrams, M. and Viskanta, R., 1974, ``The Effects of Radiative Heat Transfer Upon the Melting and Solidification of Semitransparent Crystals," Journal of Heat Transfer, Vol. 96, pp. 184-190.
Atherton, L.J., Derby, J.J., and Brown, R.A., 1987, ``Radiative Heat Exchange in Czochralski Crystal Growth," Journal of Crystal Growth, Vol. 84, pp. 57-78.
Brandon, S., and Derby, J.J., 1991, ``Internal Radiative Transport in the Vertical Bridgman Growth of Semitransparent Crystals,'' Journal of Crystal Growth, Vol. 110, pp. 481-500.
Brandon, S. and Derby, J.J., 1992, ``Heat Transfer in Vertical Bridgman Growth of Oxides: Effects of Conduction, Convection, and Internal Radiation," Journal of Crystal Growth, Vol. 121, pp. 473-494.
Cockayne, B, Chesswas, M. and Gasson, D.B., 1969, ``Facetting and Optical Perfection in Czochralski Grown Garnets and Ruby," Journal of Material Science, Vol. 4, pp. 450-456.
Crochet, M.J., Dupret, F., Ryckmans, Y., Geyling, F.T. and Monberg, E.M., 1989, ``Numerical Simulation of Crystal Growth in a Vertical Bridgman Furnace," Journal of Crystal Growth, Vol. 97, pp. 173-185.
Dupret, F., Nicodeme, P. and Ryckmans, Y., 1989, ``Numerical Method for Reducing Stress Level in GaAs Crystal," Journal of Crystal Growth, Vol. 97, pp. 162-172.
Dupret, F., Nicodeme, P., Ryckmans, Y., Wouters, P. and Crochet, M.J., 1990, ``Global Modeling of Heat Transfer in Crystal Growth Furnaces," International Journal of Heat and Mass Transfer, Vol. 33, No. 9, pp. 1849-1871.
Kvapil, Ji, Kvapil, Jo, Manek, B., Perner, B., Autrata, R. and Schauer, P., 1981, ``Czochralski Growth of YAG: Ce in a Reducing Protective Atmosphere," Journal of Crystal Growth, Vol. 52, pp. 542-545.
Matsushima, H., and Viskanta, R., 1990, ``Effects of Internal Radiative Transfer on Natural Convection and Heat Transfer in a Vertical Crystal Growth Configuration," International Journal of Heat and Mass Transfer, Vol 33, No. 9, pp. 1957-1968.
Modest, M.F., 1988, ``Radiative Shape Factors Between Differential Ring Elements on Concentric Axisymmetric Bodies," AIAA Journal of Thermophysics and Heat Transfer, Vol. 2, No. 1, pp. 86-88.
Modest, M.F., 1993, Radiative Heat Transfer, Mc-Graw Hill, New York.
Motakef, S. and Witt, A. F., 1987, ``Thermoelastic Analysis of GaAs in LEC Growth Configuration," Journal of Crystal Growth, Vol. 80, pp. 37-50.
Naraghi, M.H.N., Chung, B.T.F. and Litkouhi, B., 1988, ``A Continuous Exchange Factor Method for Radiative Exchange in Enclosures With Participating Media," Journal of Heat Transfer, Vol. 110, No. 2, pp. 456-462.
Neghabat, A. and Naraghi, M.H.N., ,1998"Analysis OF Radiative Heat Transfer in Three-Dimensional Nonhomogeneous, Nongray, and Anisotropic Scattering Media," submitted for presentation at the 11th International Heat transfer conference, Korea.
Nunes, E.M., Naraghi, M.H.N., Zhang, H. and Prasad, V., 1996, ``Combined Radiative-Convective Modeling for Material Processes: Application to Crystal Growth'', ASME Publication HTD-Vol. 323, presented at the 1996 National Heat Transfer Conference, Houston, Texas.
Nunes, E.M., and Naraghi, M.H.N., 1997, ``A Model for Radiative Heat Transfer Analysis in Arbitrarily-Shaped Axisymmetric Enclosures," ASME publication HTD-Vol. 345, presented at the 1997 National Heat Transfer Conference, Baltimore, Maryland, also Numerical Heat Transfer (in press).
Nunes, E.M., Naraghi, M.H.N., Zhang, H., and Prasad, V., 1998, ``A Volume Radiation Heat Transfer Model for Czochralski Crystal Growth Process'', Submitted for presentation at the 1998 National Heat Transfer Conference.
Nunes, E.M., Naraghi, M.H.N., Zhang, H. and Prasad, V., ``An Exchange Factor Formulation for Spectral Internal Radiation Heat Transfer Analysis of Czochralski Crystal Growth Processes," submitted for presentation at the 1998 International Crystal growth Conference to be held in Jerusalem Israel.
O'Hara, S., Tarshis, L.A. and Viskanta, R., 1968, ``Stability of the Solid/Liquid Interface of Semitransparent Materials," Journal of Crystal Growth, Vol. 3, pp. 583-593.
Rea, S.N., 1981, ``Czochralski Silicon Pull Rate Limits," Journal of Crystal Growth, Vol. 54, pp. 267-274.
Sackinger, P.A., Brown, R.A. and Derby, J.J., 1989, ``A Finite Element Method for Analysis of Fluid Flow, Heat Transfer and Free Interfaces in Czochralski Crystal Growth," International Journal for Numerical Methods in Fluids, Vol. 9, pp. 453-492.
Shapiro, A.B., 1983, ``FACET, A Radiation View Factor Computer Code for Axisymmetric, 2D Planar and 3D Geometries With Shadowing," Lawrence Livermore Laboratories, UCID-19887.
Thomas, P.D., Derby, J.J., Atherton, L.J., Brown, R.A. and Wargo, M.J., 1989, ``Dynamics of Liquid-Encapsulated Czochralski Growth of Gallium Arsenide: Comparing Model With Experiment," Journal of Crystal Growth, Vol. 96, pp. 133-152.
Tsukada, T., Kakinoki, K., Hozawa, M., and Imaishi, N., 1995, ``Effect of Internal Radiation within Crystal and Melt on Czochralski Crystal Growth of Oxide," International Journal of Heat and Mass Transfer, Vol. 38, No. 15, pp. 2707-2714.
Williams, G. and Reusser, R.E., 1983, ``Heat Transfer in Silicon Czochralski Crystal Growth," Journal of Crystal Growth, Vol. 64, pp. 448-460.
Xiao, Q. and Derby, J.J., 1993, ``The Role of Internal Radiation and Melt Convection in Czochralski Oxide Growth: Deep Interfaces, Interface Inversion, and Spiraling," Journal of Crystal Growth, Vol. 128, pp. 188-194.
Xiao, Q. and Derby, J.J., 1994, ``Heat Transfer and Interface Inversion during the Czochralski Growth of Yttrium Aluminum Garnet and Gandolinium Gallium Garnet," Journal of Crystal Growth, Vol. 139, pp. 147-157.
Yuferev, V.S. and Vasil'ev, M.G., 1987, ``Heat Transfer in Shaped Thin-Walled Semi-transparent Crystals Pulled From the Melt," Journal of Crystal Growth, Vol. 82, pp. 31-38.
Zhang, H. and Prasad, V., 1995a, Reference Manual for MASTRAPP2d, State University of New York at Stony Brook.
Zhang, H. and Prasad, V., 1995b, ``A Multizone Adaptive Process Model
for Low and High Pressure Crystal Growth," Journal of Crystal Growth, Vol.
155, pp. 47-65.