Heat and mass transfer during crystal growth

COMPUTATIONAL MATERIALS SCIENCE

ELSEVIER

Computational Materials Science 10 (1998) 127-133

Heat and mass transfer during crystal growth Koichi Kakimoto *, Hiroyuki Ozoe Institute of Advanced Material Study, Kyushu University, 6-1 Kasyga-Koen,

Kasuga, 816 Japan

Abstract Quality of semiconductor and oxide crystals which are grown from the melts plays an important role for electronic and/or optical devices. The crystal quality is significantly affected by the heat and mass transfer in the melts during crystal growth in a growth furnace such as Czochralski or horizontal Bridgman methods. This paper reviews the present understanding of phenomena of the heat and mass transfer of the melts, especially instability of melt convection from the detailed numerical calculation, which helps to understand the melt convection visualized using X-ray radiography. Large scale simulation of melt convection during crystal growth is also reviewed. Characteristics of flow instabilities of melt convection with a low Prandtl number (ratio between momentum and thermal diffusivities) are also reviewed by focusing on the instabilities of baroclinic, the Rayleigh-Benard and the Marangoni-Benard, from the points of view of temperature, rotating and/or magnetic field effects during crystal growth. Oxygen concentration in grown crystals is also discussed how melt convection affects. Copyright 0 1998 Elsevier Science B.V.

1. Introduction The Czochralski (CZ) [l] and/or horizontal Bridgman (HB) [2] crystal growth technique are widely accepted for fabricating high-quality substrates for silicon (Si) VLSIs, gallium arsenide (GaAs) monolithic and integrated circuit devices, and optical devices. It is well known that the breakdown voltage of oxide layer grown on Si substrates strongly depends on the growth conditions under which the Si crystals were grown, such as the crystal pulling speed [3] or the temperature distribution in the growth furnace [4]. The melt convection should be controlled from the point of view mentioned above. However, an actual flow has been hard to be monitored for semiconductor or oxide melts on account of opaque melts. * Corresponding author. Tel.: +81 92 583 7836; fax: +81 92 583 7838; e-mail: [email protected].

So far, the distribution of impurities and point defects in grown crystals, which affects the degradation of the breakdown voltage for silicon devices, is dependent on the amplitude of temperature fluctuation at solid-liquid interface [3]. The fluctuation is mainly caused by the flow instability which contains laminar or turbulent flows. Therefore, the origin of the flow instability should be clarified, so that it can be controlled to obtain high-quality crystals. Several possibilities of the instabilities have been reported in the CZ crucible, for example Rayleigh-Benard and baroclinic instabilities by using model fluids [5]. Visualization using numerical simulation and experimental setup would be essential for correct understanding of the flow of actual molten silicon itself to confirm the origin of the instabilities. The purpose of the present paper is to introduce hierarchy of possible flow instabilities in an actual CZ system which were clarified by a large scale computation

0927-0256/98/$19.00 Copyright 0 1998 Elsevier Science B.V. All rights reserved PII SO927-0256(97)00090-6

K. Kakimoto.

128

of fluid flow and experimental ray radiography.

Additionally,

H. Ozoe/Compututional

observation magnetic

on melt flow and oxygen distribution

using X-

in the melt are

cretizing

the governing

Navier-Stokes,

method

was used for dis-

equations

such as continuity,

energy, and impurity

transfer

formed by using alternating directional implicit (ADI) method [6] as a matrix solver. Governing equations are expressed in Eqs. (l)-(3). When the effect of axial magnetic fields was taken into account, Lorentz force

g

+

Table I Thermophysical

properties

for numerical

simulation

Density

2520

(Kg/m’)

Heat capacity

2.39 x lo6

(J/m3 K)

Dynamic

viscosity

Thermal

expansion

7 X 10-4

(Kg/ms) coefficient

(K-’

1

Electrical

conductivity

1.4 x lo-” 4.14 x 10” 168.5 45 0.3 12.3 x lo5

(S/m)

gas was adopted as expressed by Eqs. (6) [7] and (7) ]81.

0 = 3.99 x 1O23exp(-2.0

x 104/T) atoms/cm3,

in Eq. (3) as

(6) 4 = h(O(melt)

vpu = 0,

- O(gas)),

(7)

(1) where h and 0 are mass transfer coefficient

= V(l-VCD) + S@ (@ = u, T. c),

-vp + F(=

(2)

pg) + f.

(3)

fL = J x B, J = a(-V’P

(4) + u x B).

(5)

concentration in the melt, respectively. h and O(gas) are imposed to 2 and Oatoms/cm.‘, respectively [8]. Thermophysical properties used in the calculation are listed in Table 1. Numerical simulation with a three-dimensional configuration with a grid of 60 x 40 x 40 in r, 0, z directions was carried out. The geometry of the present calculation

where p, T, r are density, temperature

and diffusiv-

ity for the variables such as velocity, temperature impurity

at an in-

terface between the melt and ambient gas, and oxygen

Pg+U(PV@) SfzJ=

[ 121

equa-

tions in the present calculation. Time-dependent calculation with three-dimensional geometry was per-

(f~) expressed by Eq. (4) was included an external force,

10 (1998) 127-133

Latent heat of fusion (J/m”) Melting temperature (K) Thermal conductivity (W/mK) Emissivity

2. Numerical techniques and models volume

Science

field effects

also reviewed.

The control

Muterials

concentration.

p, f and g are pressure,

and

was set identical

to the experimental

one

with a 3 in. diameter crucible and 1.5 in. crystals. The temperature

boundary

conditions

of the melt, which

an

strongly affect the flow mode [9], were set to be axi-

external force and gravitational acceleration, respectively. U, T and c are velocity, temperature and impurity concentration, respectively. S is a source term of each valuable of velocity, temperature and impurity concentration. J, 9, CTand B are electric current, scalar potential, electric conductivity and magnetic field, respectively. When oxygen transfer in silicon melt was taken into account in the calculation, the following assumption of equilibrium concentration at an interface between the melt and a crucible, and flux at the boundary melt-

symmetric to identify the origin of non-axisymmetric profile. One of the authors has succeeded in the direct observation of molten silicon convection by using X-ray radiography method [lo] including an effect of magnetic fields [ 111. Two sets of X-ray sources and X-ray cameras are set to a furnace to visualize melt flow using a special tracer [lo]. Flow velocity obtained from the visualization can be compared with numerical result which can check a validity of both experimental and numerical results.

K. Kakimoto, H. Ozoe/Computational Materials Science 10 (1998) 127-133

3. Temperature and rotation effects Flow visualization was carried out using the following condition to clarify temperature effects. Crystal and crucible rotated with angular velocities of + 1 and - 1 rpm, respectively, to study an effect of natural convection. Crystal growth was stopped to remove a modification of temperature distribution along boundaries. Figs. l(a) and (b) show plane views of calculated velocity vectors and temperature contours on 5 mm below the surface plane. If the flow mode is exactly axisymmetric azimuthal, velocity component should not be observed. However, alternative modulation of azimuthal velocity component which is almost one

cdsec

Fig. 1. Plane views of calculated velocity vectors (a) and temperature contours (b) on the plane of 5 mm below the surface.

129

order of magnitude smaller than radial velocity existed. The lines indicated by A and B in Fig. l(a) are corresponding to the lines indicated by A and B in Fig. l(b) [12]. It can be easily recognized that fluid is ascending from inside of the melt and spreading out in the azimuthal direction along line A at the top of the melt, although it is descending into inside of the melt along line B. Lines A-C and B-C correspond to positions in which flow is ascending and descending along the lines, respectively. Additionally, temperature in the ascending part is higher than that in the descending part. These results lead us to understand that cell structure is originated by the Rayleigh and/or thermocapillary Benard instabilities [ 121. Almost the same velocity profile can be observed by experiments using X-ray radiography [lo]. From the above results, we are able to conclude that the spoke pattern is formed in silicon melt in a CZ crucible. Additionally, we can conclude that the pattern penetrates into the bottom of the silicon melt although it has been thought to be terminated just beneath the surface for the case of oxide melts [ 131. The reason for the increase of penetration depth may be attributed to the large flow velocity which transfer momentum of the melt readily to the bottom. Since silicon melt has small viscosity, the Rayleigh number (Ra = g/l ATL3/ffv which contains viscosity becomes large, where /3, AT, L, cx and v are volume expansion coefficient, temperature difference between crystal and crucible wall, crucible radius, thermal diffusivity and kinematic viscosity, respectively. Consequently, the flow becomes inertial flow dominant. This means that the large momentum of fluid transfers from surface to the bottom readily. Therefore, the pattern reaches from the surface to the bottom. The origin of the pattern is thought to be both the Benard [ 121and/or thermocapillary [12] instabilities, although it is difficult to identify which instability is dominant because reported values of temperature gradient of surface tension which originates surface-tension-driven flow have some sort of ambiguity. Therefore, it is important to obtain reliable value of the temperature gradient of surface tension which produces thermocapillary flow. We previously reported that the flow becomes axisymmetric at a AT (= 55 K) for the small crucible

K. Kakimoto, H. Ozoe/Computational

130

case: H = 3.84cm,

wC = -1 rpm, ws = 1 rpm, R, =

Materials Science 10 (1998) 127-133

respectively. Wavy structure with wave number of 2 in

3.75 m, and R, = 1.75 cm [ 141, where H, wc, us, R,,

azimuthal direction can be recognized

and R, are melt height, crucible

be observed from experiments

and crystal rotation

which can also

of X-ray radiography.

rates, crucible radius and crystal radius, respectively.

It is also clarified that the vortices are penetrated

In the second calculation,

the bottom.

both crucible and crystal ro-

tation rates were set to oC = -4rpm,

ws = -4rpm

Benard convention

into

was found in the melt which is

to study rotational effect. The calculated velocity vec-

attributed to temperature

tors from rotational

crucible and crystal rotation rates were small. Vortice’s

viewpoint

with the same angular

distribution

which was observed

in the melt when

velocity as crucible rotation rate and from stationary

structure

from a rotating

view

point of view are shown in Figs. 2(a) and (b) [14],

was observed when crucible rotation rate was large. When crystals were grown under axisymmetric flow, small inhomogeneity of oxygen was found in grown crystals. However, large fluctuation of oxygen was observed in grown crystals when asymmetric flows with vortices were formed in the melt [15]. Therefore, it is important to control melt flow to obtain homogeneous

impurity distribution

in grown crystals.

(4 4. Magnetic field’s effect Fig. city

3 shows

and

relationship

strength

cles and broken

lines represent

I

0

Fig. 2. Calculated velocity vectors from with the same angular velocity as crucible vectors trom stationary viewpomt (b).

rotational viewpoint rotation rate (a) and

of magnetic

between fields.

flow Solid

experimental

1

cirresult

I

0.02 Magnetic flux density : B (T)

Fig. 3. Flow velocity of silicon as a function strength. Solid and broken lines are calculated results, respectively.

velo-

0.04

of magnetic field and experimental

K. Kakimoto, H. Ozoe/Computational

and numerically calculated results by using threedimensional calculation [ 161. The calculated results agree well with experimental one except at low magnetic fields. The discrepancy in low magnetic fields may be due to the following reason. We obtained a steady solution from numerical calculation even if we carried out transient calculation, however, the actual flow was time-dependent. This means the present calculation offers too stabilized solution. Therefore, it is necessary to use an algorithm with higher order accuracy for non-linear term in Navier-Stokes equation. Trend of the reduction of flow velocity sometimes has been explained by a concept of “magnetic viscosity” [ 171, however, the concept does not directly express fluid motion. Fluid motion of molten silicon is able to express by “magnetic number” [17] contrary to the magnetic viscosity [ 171. The concept is based on the inertia! flow regime, therefore, thickness of velocity boundary layer near a crucible wall is almost constant with and without magnetic fields. Temperature at the bottom of the crucible was set at two different values, 1412°C (melting point) which is identical to types A in Fig. 4(a), and 1430°C which is identical to type B in Fig. 4(b) to clarify how Benard

Materials Science 10 (1998) 127-133

131

convection occurs under the two different temperature boundary conditions of types A and B [ 16]. These two types of temperature distribution could be obtained by modification of heater system shown in Fig. 4. An axisymmetric flow pattern was obtained from the numerical simulation with type A heating system within the magnetic field from 0 to 0.3T. However, a non-axisymmetric flow pattern was observed in the simulation with type B heating system as shown in Figs. 5(a) and (b), which indicate profiles of the velocity and temperature distribution at the top of the melt under a magnetic field of 0.1 T. This pattern can be found in the range of magnetic fields above 0.1 T.

8

A

(4

I

(V

(4

I

TOP

4 BOTTOM

Fig. 4. Schematic diagram of two kinds of heating system: (a) type A; (b) type B.

Fig. 5. Calculated velocity vectors (a), temperature distribution (b) and velocity vectors in a z-6 plane (c).

132

K. Kakimoto,

When magnetic velocity

becomes

Materials

Science

10 (1998) 127-133

field becomes larger than 0.1 T, flow small, therefore,

file becomes almost axisymmetric. non-axisymmetric

H. Ozoe/Computational

temperature To understand

prothe

structure in more detail, the veloc-

ity profile in z-0 plane is shown in Fig. 5(c) at a radius of 1 cm under a magnetic structure

field of 0.1 T. A cell

similar to Benard cells can be recognized

[16]. This cell structure

can be observed

only with

the type B heating system, and it was not observed with the type A system. This suggests that the origin of the cell structure was Benard instability. Since vertical magnetic fields were applied to the melt, radial flow was suppressed due to the Lorentz force [ 181. Consequently, the temperature gradient in radial and vertical directions of the melt increases and the system becomes unstable because of the cold melt sitting

u

10170

on top of the hot melt. Formation of the Benard cells relaxes the unstable temperature distribution, and the

mental

0.2

Magnetic flux density:

system becomes more stable. Calculated oxygen concentrations as a function of applied magnetic fields in the center of the crystal grown from the melt with a height of 3.75 cm are indicated

0.1

I

0.3 B(T)

I

(4 I

in Fig. 6(a) by a solid line. The experiresults are indicated

by a broken

line. The

absolute values of the numerical results are slightly different from those of the experimental results. This discrepancy may be due to unknown parameters such as the segregation coefficient, evaporation rate at the free surface, and the dissolution rate at melt-crucible interface used in the present numerical simulation. An anomaly on the oxygen concentration for both the experimental and numerical results was observed at 0.1 T. The anomaly was not observed for type B heating system as shown in Fig. 6(b). The solid line indicates calculated results and the broken line indicates experimental results. The anomaly can be attributed to the formation of Benard cells, since the strength of the magnetic field in which the anomaly was observed was almost identical to that in the formation of Benard cells. Since melt with high oxygen concentration transferred from the bottom to a solid-liquid interface, high oxygen concentration could be observed at a condition of anomaly. Oxygen concentration decreased again above magnetic fields because the oxygen diffuses to the top of the melt and evaporates into the gas phase.

1017

0

0.1

0.2

0.3

Magnetic flux density: B(T)

(b) Fig. 6. Calculated oxygen concentration as a function of magnetic field for type A (a) and type B (b).

Magnetic fields have a lot of possibilities to control melt flow, however, flow structure can modify oxygen concentration grown crystals. Therefore, it is important to clarify melt flow more quantitatively using numerical simulation.

K. Kakimoto, H. Ozoe/Computational

5. Conclusion

Large scale computation with time-dependent and three-dimensional including global model [ 19-211 is going to become a powerful tool to investigate fluid flow during crystal growth. The computation can offer more accurate result if it can use more appropriate algorithm for an actual flow of the melts of semiturbulent flow. Therefore, new algorithms are expected which can express flow structure between laminar and turbulent flow.

Acknowledgements The authors would like to acknowledge Mr. Minoru Eguchi and Dr. Yi Kung-Woo of Fundamental Research Laboratories of NEC Corporation for their skillful experiment and fruitful discussion. A part of this work was conducted as JSPS Research for the Future Program in the Area of Atomic-Scale Surface and Interface Dynamics.

References [l] J. Czochralski, Z. Phys. Chem. 92 (1917) 219. [2] W. Dietze, W. Keller and A. Muhlbauer, Crystals - Growth, Properties and Applications, ed. J. Grabmaier, Vol. 5 (Springer, Berlin, 1982).

Materials Science IO (1998) 127-133

133

[31 H. Oya, Y. Horioka, Y. Furukawa and T. Shingyoji, in: Extended Abstracts of the 37th Spring Meeting of The Japan Society of Applied Physics and Related Society (1990). r41 E. Domberger and W.V. Ammon, J. Electrochem. Sot. 143 (1996) 1636. VI J.R. Ristorcelli and J.L. Lumley, J. Crystal Growth 116 (1992) 447. [61 D.W. Peaceman and H.H. Rachford, J. Sot. Indust. Appl. Math. 3 (1955) 28. 171 H. Hirata and K. Hoshikawa, J. Crystal Growth 96 (1989) 747. PI K. Kakimoto, Y.-W. Yi and M. Eguchi, J. Crystal Growth 163 (1995) 238. [91 M. Mihelcic, C. Wingerath and Chr. Pirpon, J. Crystal Growth 69 (1984) 473. [lOI K. Kakimoto, M. Eguchi, H. Watanabe and T. Hibiya, J. Crystal Growth 88 (1988) 365. [Ill K. Kakimoto and K.-W. Yi, Physica B 216 (1996) 406. [I21 K.-W. K, K. Kakimoto, M. Eguchi, M. Watanabe, T. Shyo and T. Hibiya, J. Crystal Growth 144 (1994) 20. [I31 A.D.W. Jones, J. Crystal Growth 65 (1983) 124. [I41 K.-W. Yi, V.B. Booker, M. Eguchi, T. Shyo and K. Kakimoto, J. Crystal Growth 156 (1995) 383. [151 M. Watanabe, M. Eguchi, K. Kakimoto, H. Ono, S. Kimura and T. Hibiya, J. Crystal Growth 151 (1995) 285. [I61 K. Kakimoto, K.-W. Yi and M. Watanabe, J. Crystal Growth 163 (1996) 238. [I71 K.-W. Yi, M. Watanabe, M. Eguchi, K. Kakimoto and T. Hibiya, Jpn. J. Appl. Phys. 33 (1994) L487. [181 H.P. Utech and MC. Fleming, J. Appl. Phys. 37 (1966) 2021. [I91 R.A. Brown, T. Kinney, P. Sackinger and D. Bomside, J. Crystal Growth 97 (1989) 99. WI F. Dupret, P. Nicodeme, Y. Ryckmans, P Wouters and M.J. Crochet, Int. J. Heat and Mass Transfer 33 (1990) 1849. [211 A. Virzi and M. Porrini, Mat. Sci. Eng. B 17 (1993) 196.