Dopant transport during semiconductor crystal growth with magnetically damped buoyant convection

,. . . . . . . .

ELSEVIER

CRYSTAL QROWTH

Journal of Crystal Growth 172 (1997) 124-135

Dopant transport during semiconductor crystal growth with magnetically damped buoyant convection N. Ma, J.S. Walker * Department of Mechanical and Industrial Engineering, Unil,ersity of Illinois, Urbana, Illinois 61801, USA

Received I 1 June 1996;accepted 5 August 1996

Abstract

This paper presents an asymptotic model for the unsteady transport of a dopant during the growth of a semiconductor crystal from a melt with an externally applied magnetic field. The melt is divided into (1) mass-diffusion boundary layers where convective and diffusive mass transfer are comparable, and (2) a core region where diffusion is negligible, so that the concentration of each fluid particle is constant. A Lagrangian description of motion is used to track each fluid particle during its transits across the core between diffusion layers. The dopant distribution in each layer depends on the concentrations of all fluid particles which are entering this layer. The dopant distribution is very non-uniform throughout the melt and is far from the instantaneous steady state at every stage during crystal growth. The transient asymptotic model predicts the dopant distribution in the entire crystal.

I. I n t r o d u c t i o n

Turbulent and unsteady laminar melt motions produce undesirable spatial oscillations of the dopant concentration (striations) in a crystal grown from a melt. Since most molten semiconductors are good electrical conductors, an externally applied magnetic field can eliminate turbulence and all other unsteadiness in the melt motion, thus eliminating the striations in the crystal. The first papers to demonstrate the benefits of applying magnetic fields during semiconductor crystal growth appeared in 1966 [1,2]. Unfortunately the elimination of turbulent mixing and a moderate electromagnetic (EM) damping of the residual steady laminar melt motion may lead to a large variation of the crystal's dopant concentration * Corresponding author. Fax: + 1 217 244 6534.

in the direction perpendicular to the direction of growth (lateral or radial macrosegregation) [3]. If the EM damping is extremely strong, then the melt motion has no effect on the dopant distribution in the crystal, and this diffusion-controlled mass transport may produce a dopant distribution in the crystal which is both laterally and axially uniform [3]. In order to achieve diffusion-controlled dopant transport, the mass transport Prclet number, Pe,~ = U b L / D , must be small, where U b is the characteristic velocity for the electromagnetically damped melt motion and varies as some negative power of the magnetic flux density, while L and D are the characteristic length for the melt and the diffusion coefficient for the dopant in the molten semiconductor. Since typical values of D are 1 to 2 × 10 -8 mZ/s, Ub must be incredibly small for diffusion-controlled dopant transport, and the required magnetic flux

0022-0248/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S0022-0248(96)00720-8

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

density may not be possible or practical. The reasonable limit of current superconducting magnets is about six tesla (T). The cost of a magnet is roughly proportional to the magnetic flux density times the volume, and large volumes are required for crystal growth because the entire furnace must fit inside the magnet. Therefore it is more practical to use magnetic fields which are strong enough to eliminate turbulence and other flow unsteadiness, but which only moderately damp the residual laminar melt motion, i.e., for which Pem is still large. With such fields, the objective is to tailor the melt motion in order to achieve both lateral and axial dopant uniformity in the crystal. This is the objective of the non-uniform, axisymmetric or " c u s p " fields currently being used for the Czochralski growth of silicon crystals [4]. With this tailored magnetic field, the axial field at the crucible bottom and the radial field at the vertical crucible wall strongly suppress the buoyant convection which carries oxygen from the crucible to the melt, the radial field at the free surface weakly suppresses the thermocapillary convection needed to insure that most of the oxygen evaporates, and the radial field at the crystal-melt interface weakly suppresses the radial flow due to the crystal rotation needed to produce a radially uniform oxygen distribution in the crystal. At each stage during the growth of a crystal by any process, there are infinitely many different ways to tailor the strength and configuration of the externally applied magnetic field, the rotation rates of the crystal and crucible or ampoule or feed rod, the distribution of the heat flux into the melt, the radiative and conductive heat losses from the melt, etc., so that process optimization through trial-and-error experimental crystal growth is not practical. Models which accurately predict the dopant distribution in an entire crystal for any combination of process variables are needed to facilitate process optimization. Most previous treatments of dopant transport have involved the numerical solution of the full mass transport equation over the entire melt [3,5,6]. With a numerical treatment of the full mass-transport equation for an unsteady dopant concentration in the melt, the time step must be extremely small and the spatial grid must be extremely fine because the small value of D leads to very thin mass-diffusion boundary layers [5]. Therefore hundreds of hours of super-

125

computer time are needed to simulate five to ten minutes of transient dopant transport [5]. It is clearly impractical to extend this transient numerical approach to the period of twelve hours to several days needed to grow a single crystal. One alternative is to treat the dopant distribution in the melt as an instantaneous steady state at each stage during the growth of a crystal [3,6]. The assumption that the dopant distribution reaches an instantaneous steady state is probably valid for no magnetic field and for the extremely strong magnetic fields producing diffusion-controlled mass transport, but it is not valid for the moderate field strengths which eliminate turbulence and unsteadiness in the flow, but which lead to moderate or large values of Pe m. Without a magnetic field, turbulent mixing greatly accelerates the approach to an instantaneous steady state. For diffusion-controlled mass transport, none of the rejected dopant escapes from the boundary layer adjacent to the crystal-melt interface, and this layer quickly approaches a steady state. With moderate magnetic field strengths, the residual laminar melt motion convects the rejected dopant out of the crystal-melt boundary layer, producing very non-uniform and intrinsically unsteady dopant distributions throughout the melt. This paper presents an asymptotic treatment of the unsteady dopant transport for the entire period needed to grow a crystal. We use a two-dimensional model problem with a simple melt motion in order to illustrate the asymptotic method, but the extension to the three-dimensional melt motions in actual crystal growth processes would be straightforward, and the required computational resources would still be quite modest. Series and Hurle [4] review the literature on semiconductor crystal growth in magnetic fields, with an emphasis on the mass transport of dopants and contaminants. Our Pe m is the same as the parameter ( S c G r / H a 2) used by other researchers [7], where Sc, Gr and Ha are the Schmidt, Grashoff and Hartmann numbers, respectively, 2. Melt motion

This paper treats the unsteady, two-dimensional mass transport of a dopant in a solidifying, electrically conducting semiconductor melt in a horizontal,

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

126

~3

INTERFACE x 2

._......._.-~

MELT

-.--------.--.~ _ _

CRYSTAL

[ g

< a

Fig. 1. Two-dimensional model problem with a unilbrm, steady, vertical magnetic field B0.f and with coordinates normalized by half the distance between the top and bottom walls. rectangular container with thermally insulated top and bottom walls and with an externally applied, uniform, vertical magnetic field, B0.~. Here B 0 is the magnetic flux density, while .f, ,f and ~ are unit vectors for the Cartesian coordinate system. Our dimensionless problem is sketched in Fig. 1. The coordinates and lengths are normalized by half the distance between the top and bottom walls L, and a is the dimensionless length of the container. The melt velocity is normalized by the characteristic velocity for magnetically suppressed buoyant convection [8], pg /3qL ub ' (1) where p is the melt density at the solidification temperature ]C,, g = 9.81 m / s 2, and q is the constant, uniform density of the heat flux into the right end of the container, while /3, o- and K are the volumetric expansion coefficient, electrical conductivity and thermal conductivity of the melt, respectively. The density of the heat flux removed at the left end of the container is adjusted so that the crystal-melt interface moves at a constant velocity, Ug = (.DUb. With time t normalized by L / U b, the dimensionless time to grow the entire crystal is O/O).

The electric current in the melt produces an "induced" magnetic field, which is superimposed on the "applied" field produced by the external magnet. The characteristic ratio of the induced to applied magnetic field strengths is the magnetic Reynolds number R m = ],~pO'UbL ,

(2)

where /% is magnetic permeability of the melt. For all crystal growth processes, R m << 1, and the induced magnetic field is negligible. The characteristic ratio of convective heat transfer to conduction is the thermal P6clet number, Pe t = pChUbL/K, where c h is the specific heat of the melt. For a sufficiently strong magnetic field, Pe t << 1, convective heat transfer is negligible, and the crystal-melt interface is a vertical plane, as shown in Fig. I. In general, Ug_> Ha 3/2 [91. In our two-dimensional model problem, nothing varies in the z-direction, and there can be a uniform electric field in the z-direction E_. In an actual horizontal Bridgman process, there are electrically insulating walls at, say, z = _ d , which block any

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

127

Harl h

h

,,/ Y

I I

I I

T

x

~

H a "1/2 I Ha'l/2

3,

.

I

, I

I x=~t

--

--

--I

x=a

' / L__,=.,

~ Ha-1

h

Fig. 2. Flow subregions of the melt for Ha >> 1: c = inviscid core, h = Hartmann layers adjacent to the top and bottom, and p = parallel layers adjacent to the crystal-melt interface and right end of container, which are parallel to the magnetic field.

net electric current in the z-direction. For the present recirculating flow, E. = 0 for zero net current. O h m ' s law gives j = u£ for the electric current density normalized by o-U b B0, where v = u2 + @ is the dimensionless melt velocity. With the Boussinesq approximation, the continuity and inertialess momentum equations are V "t, = 0,

zc

~0(s~,y) = ½ ( 1 - y 2 ) -

Ha-2V2v,

,o.

(Sa)

0q, Vp, = - Ha'/2 07 ( ~:, y ) ,

Z,cos[(n+½)Try]

X exp( - a . sc) [cos( a. sc) + sin( a . ~ ) ] ,

(5c)

(4b)

where p is the deviation of the pressure from the hydrostatic pressure for the uniform density p, normalized by p g ~ S q L 2 / K , while - u 2 is the EM body force due to the electric current in the z-direction and T in the buoyancy term is given by Eq. (3). For Ha >> 1, viscous effects are confined to thin boundary layers and the melt can be divided into the subregions shown in Fig. 2: (c) is the inviscid core, (h) are Hartmann layers with A y = O(Ha -t ) at y = ___1, and (p) are parallel layers with A x = O ( H a - J/2 ) at x = wt and at x = a. In the inviscid core region, Pc = ( x - o)t)y, u c = - y and c c = 0 to all orders in Ha. Since the flow circuit must be completed through the two parallel layers, u = O(1) and v = O(Ha 1/2) inside these layers. Inside the left parallel layer adjacent to the crystal-melt interface, the velocity measured in a frame of reference moving with the interface is -

• n=0

(4a)

0 = - V p - u.f + T~ +

.p, =

neglecting O ( H a - L) and O ( H a - i/2) terms, respectively. Here, ~:= H a l / 2 ( x - ~ot), and the separationof-variables solution for the stream function $ is [101

(5b)

An = ½['rr(2n + 1)] ,/2, A,+=16(-1)n[~r(2n

+ l)]

(5d) 3

(5e)

Inside the right parallel layer, the velocity measured in the fixed frame is

00 Upr = ~ - 7 ( ~ ' , y ) ,

(6a)

Vor = Hal/2 0q* ( ~ ' , y ) , 0~'"

(6b)

where ~"= H a l / 2 ( a - x). The Hartmann layers have a simple exponential structure which satisfies the no-slip condition u = 0 at the top or bottom and which matches the values of u in the core or in either parallel layer. Our two-dimensional model problem involves (1) a simple linear variation of the horizontal velocity in the core, with flows to the left and right for y > 0 and y < 0, respectively, and (2) O(Ha 1/2) vertically downward and upward velocities inside the parallel layers adjacent to the crystal-melt interface and right end of the container, respectively. There are two

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

128

differences between our model problem and an actual horizontal Bridgman process with a strong magnetic field [11]. First, the heat flux into the melt in an actual process occurs across the top and bottom boundaries near the crystal-melt interface, so that (1) the melt near the hot end of the actual Bridgman boat is isothermal, and (2) the flow circuit consists of an O(Ha ~/2) vertically downward velocity inside the parallel layer adjacent to the crystal-melt interface and streamlines which close in the core near the interface. Second, if the actual Bridgman boat has vertical, electrically insulated sides at z = _+d, then there are O(Ha ~/2) velocities inside the parallel layers adjacent to the sides, and the flow is intrinsically three-dimensional. The purpose of the present paper is to illustrate the asymptotic treatment of the unsteady dopant transport during the entire period of time required to grow a crystal. The simplicity of our core solution facilitates the Lagrangian tracking of melt particles, so that we can illustrate the characteristics of the asymptotic treatment in a simple context. The extension to a realistic melt motion would simply involve the addition of an elaborate bookkeeping process to track every melt particle from the start of crystal growth until its solidification or the end of crystal growth.

at t = 0. Once crystal growth starts, rejected dopant elevates the concentration near the growth interface, and the melt motion convects the elevated concentration into the rest of the melt. The dimensionless equation governing this mass transport is

OC - - + v ' VC = PemJV2C, at

where Pe m = U b L/D is the mass P6clet number and D is the diffusion coefficient. Since dopant concentrations are generally small, the dilute approximation is appropriate. We have implicitly assumed that the dopant density has no effect on the buoyant convection. As an example, we use the properties of molten indium-phosphide, L = 1.7 cm, a temperature gradient of 19 K / c m in the melt and D = 2 X 1 0 ~ m2/s, so that Ha = 497B 0 and Pem = 854.7Bo 2. For B 0 = 4 T, Ha j/2 = 44.6 and Pe m = 53.4 so that it is appropriate to take Pem = y H a j/2,

Before solidification begins, the dopant concentration is uniform, and this initial value is used to normalize the concentration C, so that C(x, y,t)= 1

Y

lx

I

T

--I

s

,

a -IN

I

d ~

c

S

--

I

.I P a-l/4

1

d

i" h /

Ha-1/z

I

Ha'l/2 I~'~"

---4

Haq/2

__l-- - -

I

I

Haq/21 --

s

I

"t . . . .

>

Mass transport

[

,

p

x=oll

(8)

where y is an O(1) parameter. For this case, the mass-diffusion boundary layers adjacent to the crystal-melt interface and to the right end of the container have the same thickness as the parallel layers in the melt-motion solution, namely O(Pem~) = O(Ha-J/2), because the layer velocities which are perpendicular and parallel to the boundary are O(1) and O(Ha~/2), respectively. The subregions of the melt for the asymptotic treatment of the dopant transport with Pem = y H a ~/2 are shown in Fig. 3. In addition to the melt-motion subregions in Fig. 2,

3. Asymptotic solution for mass transport

Ha-it'2

(7)

~---

I

s

r--

x=a

J Ha-1/2

-----4----

y=-I

~Ha "1

Fig. 3. subregions of the melt for Pem = O(Ha I/2) and Ha >> 1. In addition to the flow subregions in Fig. 2, there are lbur square regions denoted by s and two diffusion layers denoted by d.

N. Ma, J.S. Walker/Journal of Cr3.,stal Growth 172 (1997) 124-135

there are diffusion layers with an O(Peml/2) = O(Ha -1/4) thickness between the core and its top and bottom Hartmann layers, and there are square regions with Ax = A y = O(Pem l) = O(Ha -1/2) between each parallel layer and its top and bottom Hartmann layers. The diffusion layers near y = + 1 are much thicker than the parallel layers because their velocities perpendicular and parallel to these walls are zero and O(1), respectively. In the treatment of the dopant transport in the same two-dimensional problem, but without a magnetic field, Garandet [12] found regions which are similar to our square regions and which had the dimensions A x = A y = (GrSc) -2/9. Since our Pe m = G r S c / H a 2, Eq. (8) implies that GrSc = y H a s/2, so that our square region dimensions correspond to Ha ~/2 = (GrSc) 1/5, which are slightly larger than those found by Garandet [12]. The slight difference arises because the melt-motion solution used by Garandet [12] has no thin viscous boundary layers and satisfies the no-slip condition at every boundary, while our Hartmann layers are much thinner than the square regions, so that u and u are O(1) in our square regions and the thinner Hartmann layers satisfy the no-slip conditions at y = _+ 1. For our indium-phosphide example with a typical growth rate Ug = 3 Ixm/s, the Prclet number of solidification Pe~ = UgL/D--2.55, so that we assume that Pe~ = O(1). Before we present the solution for Pe m = y H a 1/2, we discuss all seven different mass transport solutions for different relationships between Pe m and Ha, for Ha >> 1 and for Pe s = O(1). The seven different relationships are: (1) Pe m << 1, (2) Pem = O(1), (3) 1<< Pem << Ha 1/2, (4) Pe m = O(Hal/2), (5) Ha l / 2 < < P e m < < H a 2, (6) Pe m = O(Ha2), and (7) Pe m >> Ha 2. For our indium-phosphide example, these cases correspond to: (1) B 0 >> 29.2 T, (2) B 0 = 0(29.2 T), (3) 4.30 T << B 0 << 29.2 T, (4) B 0 = 0(4.30 T), (5) 0.243 T << B 0 << 4.30 T, (6) B0 = 0(0.243 T), and (7) B 0 << 0.243 T, respectively. The only case for which the melt motion has no effect on the dopant distribution in the crystal is case (1) with Pem << 1. This case can never be achieved for crystal growth on earth, but can be achieved with magnetic damping during crystal growth in space because the magnitude of the residual accelerations in an earth-orbiting vehicle is much smaller than that of terrestrial gravity [13]. For case

129

(3) with 1 << Pe m << Ha 1/2, the parallel layers in Fig. 3 are split into thinner viscous layers with O(Ha -1/2) thickness and thicker mass-diffusion layers with O(Pem ~) thickness. Correspondingly, the dimensions of the square regions in Fig. 3 are larger, namely A x = O(Pe m1) and A y = O(Pe m1). For case (5) with Ha 1/2 << Pem << Ha 2, the parallel layers are split into thicker viscous layers with O(Ha -1/2) thickness and thinner mass-diffusion layers with O((HaPem)- 1/3) thickness, while the square regions have shrunk to Ax=Ay=O((HaPem)-l/3). The only difference between case (5) and case (6) with Pe m = O(Ha 2) is that the square regions have shrunk to A x = A y = O ( H a - l ) for case (6), so that the velocity in each square region is no longer given by a Taylor series expansion of the viscous parallel-layer velocity near ~ = 0 or ff-- 0 and y = + 1. We used two procedures in order to insure that our solution for Pe m = y H a 1/2 gives accurate results for cases (3), (4) and (5), i.e., for small, moderate and large values of y, respectively. Therefore our model is valid for 1 << Pe m << Ha 2, or 0.243 T << B 0 << 29.2 T for our indium-phosphide example. Our first procedure is to insure that the numerical method adequately resolves both layers when each parallel layer is split into an O(Ha -1/2) viscous layer and an O(Pem I) or O((HaPem) -1/3) mass-diffusion layer. For 1 << Pem << Ha 1/2, i.e., small values of y, the parallel-layer velocity given by Eqs. (5) and (6) includes the core velocity, so that we need only insure that we truncate the parallel-layer computational domains at sufficiently large values of ~ and ~" in order to include all dopant diffusion. For Ha 1/2 << Pe m << Ha 2, i.e., large values of y, we need only insure adequate numerical accuracy near ~:= 0 and i f = 0 in order to resolve the relatively thin O((HaPem ) - 1/3) diffusion layers. Our spectral solution automatically concentrates the collocation points near ~ = 0 and i f = 0 , so that we can accurately cover the entire case (5) range with quite modest computational resources. Our second procedure is to use a composite solution for each parallel layer and its pair of square regions. If we use separate solutions for the parallel layers and square regions, we need a separate analysis for each of the three cases (3), (4) and (5) because the matching conditions between each parallel layer and its square regions are different for each case. However, if we use a c o m -

130

N. Ma, J.S. Walker/Journal of Crystal Growth 172 (1997) 124-135

posite solution which is valid for the parallel layer and its square regions, then no matching is required, and the solution applies equally well for all three cases. Since the vertical coordinate for the composite solution is y, we must insure that the numerical method adequately resolves the y derivatives in the square regions for all three cases. The high density of collocation points near y = _+1 in our spectral solution means that this resolution is again achieved with modest computational resources. Since Ha appears as a parameter in the boundary value problem for the composite parallel-layer-square-region solution, we must specify the value of Ha along with those of Y and co for each problem. Our approach is to treat a particular material, heat flux density and geometry, and to vary B 0, e.g., our indium-phosphide example, so that there is only one variable B 0, while Ha, T and co are known for each value of B o. Within the limits of our melt-motion solution, we could easily extend our mass-transport solution to include cases (6) and (7) as well. We would use a composite solution for the velocity in the parallel layers and adjacent Hartmann layers. This would simply add terms with e x p [ - H a ( 1 -T-y)] to the velocities given by Eqs. (5) and (6), so that the no-slip conditions are satisfied at y = _+1 as well as at ~ = 0 and ~ ' = 0 . The composite parallel-layersquare-region solution for the mass transport using the composite parallel-layer-Hartmann-layer solution for the velocity would be valid for all Pem >> I. However, the values of B 0 for cases (6) and (7) are generally sufficiently small that convective heat transfer and inertial effects are not negligible, so that our simple melt-motion solution is no longer valid and the crystal-melt interface is no longer planar and vertical. For Pe m >> 1, Eq. (7) indicates that the dopant concentration C of a fluid particle remains constant during each of its transits across the core, and that its C only changes as the particle passes through the left or right parallel layer. For each value of y, the core concentration C c is given by C~( x , y , t ) = F ( x + y t ) ,

(9)

where the function F is determined by the concentrations of the particles leaving the left and right parallel layers. For y > 0, the fluid particle entering the left parallel layer at x = cot at the present time t

left the right parallel layer at x = a at the previous time ~ = t - (a - cot)/y. For - co _< y < O, the core melt moves to the right, but it moves slower than the crystal-melt interface, so that it is overtaken by the left parallel layer. Therefore the core fluid here has never been in either layer, so that C c ( x , y , t ) = l,

for

-co_
(10)

for all t, and this fluid is entering both parallel layers. For - l < y _< - co, the particle which enters the right parallel layer at x = a at the present time t left the left parallel layer at x = co? at the previous time ? = (a + y t ) / ( c o + y). The concentrations convecting across the core do not satisfy the conditions 0C 0y

-0,

at

y=_+l,

(11)

so that the diffusion layers must satisfy these conditions and must match the non-zero values of O C J O y at y = _+ 1. For the bottom diffusion layer, we introduce the stretched coordinate v = H a l / 4 ( y + 1), and we use a reference frame which is moving with the fluid particle that exits from the left parallel layer at y = - 1 and at time t', so that ?

(t-x)

x=t-(l-co)?,

(l-co)'

for

t>?. (12)

The concentration in the bottom diffusion layer, neglecting O(Ha - j / z ) terms, is Cj = C ~ . ( x , - 1,t) + Ha - I / 4

×

[

×exp

cc

,/2 f, (;-,-31-/2

( , 2)1 4(t-

r)

dr ,

(13a)

where 0Cc ( x , G(T,?) = --~-v .

1,r) x= r (I

(13b) -

,o)?

As the melt near y = - l moves from the left layer to the right layer with the velocity u = l - Ha-~/4v, dopant diffuses in the positive y-direction, producing a more uniform vertical distribution in a thin region

131

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

whose O(Ha-~/4)Ay grows parabolically with distance from the left layer. This diffusion only involves a redistribution of the O(Ha -~/4) perturbation to the core concentration over an O(Ha -~/4) vertical distance, so that there is an O(Ha -~/2) error when we ignore this redistribution in the concentrations entering the right layer. The composite parallel-layer-square-region concentration is Cl(~,y,t) or Cr(~,y,t) for the left or right parallel layer and its pair of square regions. For Pem = y H a 1/2, the governing equations are 3C 1

3C 1 '~

3 2 C¿

T Upl----~ 3 +Ha-l/2t'pl -'~'-'1 oy 1= ~og +Ha-!

02 Cr 3~ 2

i _ _

3y2 '

3C,

t o y ( 1 - k ~ ) C l, =0,

at

~=0,

(15a)

at ,~ = 0,

(15b)

where k~ is the segregation coefficient. Matching the Hartmann layers gives

(14a)

3G

32 Cr + Ha-

3C~

32 C 1 0y2 ,

3C~ 3Cr 1 ")/ --Upr 3--~"~- H a - l / 2 Upr 3--;)

-

cles in the core because the core transit time is O(1), while the transit time through each parallel layer is negligible, namely O(Ha-l/2), because v is O(Ha 1/2) inside these layers. The boundary conditions at the crystal-melt interface and at the right end of the container are

(14b)

where upl(~,y), epl(~,y), Upr(~,y)and r'pr(~',y)are given by Eqs. (5) and (6). While the term Ha 132C/3y 2 is O(Ha -1) in each parallel layer, it is O(1) in each square region where ~/Oy = O(Hal/2). In order to obtain a composite solution which is equally valid for each parallel layer and for its pair of square regions, we only delete terms in Eq. (7) which are O(Ha -~/2) or smaller in both the parallel layer and the square region [14]. While we do not rescale the y coordinate for the square regions, we recognize that the y derivative of the composite solution becomes large near y = + I, and we insure that the numerical method resolves the large local derivative. While C~ and C r are functions of time, 3C/3t does not appear in Eqs. (14) because this term is O(Ha -~/2) for both the parallel layers and square regions. There is a transient in each parallel layer when the concentration of the entering fluid changes or in the left layer at the start of crystal growth, but the length of time needed to reach a steady-state concentration after each such change is very short, namely At = O(Ha-l/2). Therefore, while the coreconcentration problem is intrinsically unsteady, the composite-layer problems can be treated as instantaneous steady states for the instantaneous entering concentration. In a similar way, we must track parti-

3cr

. . . . 3y 0y

0,

at

y=+l.

(16)

The results of matching the parallel-layer and core concentrations depend on whether the fluid is entering or leaving the parallel layer. If the fluid is entering the parallel layer, the core concentration gives the value of C~ or C r, while the core concentration at each value of y is determined by the concentration leaving the other parallel layer at the appropriate previous time. If the fluid is leaving the parallel layer, then the parallel-layer concentration gives the value of C~, and the appropriate condition on the parallel-layer solution is that C I and C r do not grow exponentially as ~ ~ ~ or as ~"~ ~. For our numerical solution, we truncate the parallel-layer domains at ~ = ~max and ff = ~'max" There are many ways to exclude exponential growth, and we found that specifying zero gradient at the truncated boundary worked well. Therefore, the conditions at the truncated boundary are

Cl(~max,Y,t ) =Cc(tot,y,t ), C,(~max,Y,t )=

1,

for

for

0
-w
1, (17a) (lVb)

3C!

3~ ( ~ m a x , y , t ) = 0 , 3Or

ff. . . . y , t ) = 0 ,

for

-- l < y _< -- w,

(17C)

for

0
(17d)

Cr(~max,Y,t)=1, for - t o N y < 0 , Cr(~max,Y,t) =Cc(a,y,t ), for - 1

(17e)
-to. (17f)

132

N. Ma, J.S. Walker/Journal of Cr3'stal Growth 172 (1997) 124-135

At t = 0, crystal growth starts, and C~ quickly reaches the steady state corresponding to Cc = 1 in Eq. (17a). The first particle with C > I to reach the right layer arrives at t = a, so that C r = 1 for 0 < t < a. The first particle to leave the right layer with C > 1 reaches the left layer at t = 2 a / ( t o + 1), and C~ then begins to change from the solution for C c = 1 in Eq. (17a). For the Chebyshev spectral collocation solutions for C 1 and C r, we use G a u s s Lobatto collocation points in y, so that we track the particles crossing the core at the same values of y. We solve for C 1 and Cr at the discrete times t = t,, = all + n ( l t o ) / ( t o × NT)], where NT is the large number of equal-length time intervals, and n = 0, I, 2, 3 . . . . . At each time, the values of C~. in Eqs. (17a) and (17f) are determined by (1) computing the previous time t" when the particle at each value of y left the other parallel layer, and (2) interpolating between the concentrations leaving that layer at that y at the two closest previous discrete times, i.e., at ~k and ~k+l, where ?k < ~"< ?k+l. The value of NT must be sufficiently large that the minimum transit time across the core in either direction is much longer than the time interval. Since the transit time decreases as the interface moves, this scheme always fails at some time just before the end of crystal growth at t = alto. If the time interval is too large, the truncation error associated with the linear interpolation between concentrations at ~k and ~k+ t becomes significant during the final stages of crystal growth. This truncation error produces a loss of dopant, so that the accuracy at any time is easily checked since the integral of the concentration over both the crystal and the remaining melt should always equal 2 a . One could easily change the interpolation scheme or take smaller time intervals near the end, but we simply chose sufficiently large values of NT, so that the equal-time-interval scheme gives accurate results for more than 95% of the time to grow the entire crystal. Assuming that there is no dopant diffusion in the solid crystal, the dopant distribution in the crystal, C~(x,y), normalized by the initial uniform dopant concentration in the melt, is given by

C~( x, y) = k~Cl( ~ = O, y,t = x / t o ) .

(18)

For the present simple melt motion, the dimensionless length a of the container and of the entire

crystal is not a parameter. All transit times and hence all times at which given events occur are proportional to a, so that results presented as functions of t / a and x / a are valid for all values of a. The remaining parameters are k~, y, to and Ha.

4. Results We present results for k~ = 0.1 and for our indium-phosphide example with Ug = 3.0 I x m / s and with two values of B 0, so that y = 38.3Bo 5/2, w = 0.00298 B~, and Ha = 497B 0. W e set a = 1, so our x and t represent x / a and t / a . For B o = 4 T, y = 1.20, 09 = 0.0477, Ha = 1988, and the time to grow the crystal a / t o = 20.97. The constant-concentration lines in the left layer for 0 < t _< 2 a / ( 1 + to) = 1.91 are presented in Fig. 4. Fluid enters with C = I for - 0 . 0 4 7 7 < y < 1 . The fluid which solidifies at y = 1, enters at y = 1 and travels a short distance across the layer, accumulating relatively little excess dopant. The fluid which solidifies at y = - 1, enters the layer at y = 0.951, so that it must traverse the full height of the layer in the dopant-rich region near the interface, thus accumulating nearly five times as much excess dopant as the fluid at y = 1. Fluid entering the layer near y = 0.9, reenters the core near y = - 1, so that it has traveled a long vertical distance relatively close to the interface, so that it has accumulated nearly as much excess dopant as the fluid solidifying at y = - 1 , and this large dopant concentration will be carried rapidly across the core to enter the right layer at y = - 1. On the other hand, a fluid particle entering the left layer at y = 0.2, travels only a short vertical distance far from the interface, and reenters the core at y = - 0 . 2 9 5 with relatively little excess dopant. What little excess dopant it carries moves slowly across the core, so that C o ( a , - 0 . 2 9 5 , t ) = I long after C o ( a , - 1 , t ) has grown to very large values because of convection from the left layer by the high velocity near y = - 1. The fight-layer concentration when nearly 40% of the crystal is grown in shown in Fig. 5. The entering fluid has a very non-uniform concentration (1) because the concentration leaving the left layer is non-uniform, and (2) because all concentrations in the left layer increase after t = 1.91, and the fluid entering the right layer near y = - 1 or

N. Ma, J.S. Walker~Journal of Crystal Growth 172 (1997) 124-135

y

1.0

1.0 -~

0.5

0.5

0.0

y

1.05

133

0.0

~ 1 . 1 0 -0.5

-0.5 1.15

42

1.26

-I.0

i

~

0

. . . .

I

i

i

i

4

2

Fig. 4. Left layer's concentration

-ro

I

'I

6

Cl(~:,y) for

'

I

0

B0=4

for

T

4.7

'

2

Fig. 6. Left layer's concentration

I

'

4

I 6

Cl( ~: v ) f o r B o = 4 T at t = 18.7.

0_
near y = - w exited from the left layer a short or long time before, respectively. For y = 1.20, diffusion dominates over convection in the right layer so that the exiting fluid has a relatively uniform concentration. The concentration in the left layer when nearly 90% of the crystal is grown is shown in Fig. 6. While the concentration leaving the right layer at each time is relatively uniform, these concentrations are increasing with time. The fluid entering the left layer near y = 1 or near y = 0 exited from the right layer recently when the concentration was high or

1.0 -

long ago when the concentration was low, respectively. The non-uniform entering concentration cancels some of the lateral segregation evident in Fig. 4 for uniform entering concentration. The fluid which solidifies at y = 1 acquires very little additional dopant as it crosses the top of the left layer, but it had a large concentration when it entered. The fluid solidifying at y = - 0 . 5 accumulates considerable additional dopant during its long, vertical trip near the interface, but it entered with a lower concentration. Therefore there is lateral uniformity over a large part of the crystal, as shown in Fig. 7, which corresponds to Fig. 6. For - 0.5 < y < 1 in Fig. 7, there is only a 2.3% difference between the maximum and minimum values of C~. The concentration

0.5 = 0.56 0.54 y

0.0= 1.25

0.52

1.40

0.5 C

1.50

-0.5

0.48

1.60 0.46

-1.0

I 7,5

'

I 5.0

I

I

I 2,5

0.44

I 0.0

0.42 -I

Fig. 5. Right layer's concentration t = 8.09.

Cr(~,y)

for

B0=4

T

at

I

I

I

-0.5

0

0.5

Y

Fig. 7. Crystal concentration C~(0.90,y) versus y for

Bo = 4 T.

N. Ma, J.S. Walker~Journal of Co,stal Growth 172 (1997) 124-135

134

near y = - 1 in the crystal is always larger because there is a stagnation point in the flow at y = - 1 and at a small value of ~:, so that the nearly stagnant fluid near y = - 1 and sc = 0 always accumulates a large amount of dopant before solidifying. In Fig. 7, (7., at y = 1 is 25% larger than the minimum C~ = 0.433 at y = 1. The dopant distribution in the crystal, C~(x,y), for B 0 = 4 T is presented in Fig. 8. All constant-concentration lines are horizontal for 0 < x < 0.091, corresponding to the left-layer concentration in Fig. 4 for 0 < t < 1.91, and evident in the 6"., = 0.12 line in Fig. 8. The effect of dopant return to the left layer is evident in the axial increase of C~, and the effect of the concentration of the returning dopant near y = +1 is evident in the lateral uniformity over much of the crystal for x > 0. I. The lateral average of the concentration in Fig. 8 is very close to the classical well-mixed values, C~ = k~(1 - x ) (1-kO, even though the concentration in the melt is extremely non-uniform at every time. For B 0 = 0 . 5 T, "7=217, ~o=0.000745, H a = 249, and the time to grow the crystal is a / w = 1342. Of course the dimensional time to grow the crystal is the same for all values of B 0 since Ug is the same, but our characteristic time, L / U b, is proportional to Boa. With the larger value of y, convection dominates over diffusion in both layers. In the right layer, the constant-concentration lines very nearly coincide with the streamlines, so that the concentration of a particle entering the left layer is nearly equal to its value

1.0--

.) 0/00)0

1.0--

0.5,

0,18

y

-1.0

0.0-

~

' 0,00

I 0.25

0.25 0.45

1 0,50

0.75

x

Fig. 9. Crystal composition C~(x,y) for B0 = 0.5 T.

when this particle left this layer. Of course, the delay times for the return of rejected dopant are short and long for particles near y = + 1 and near y = 0, respectively. For y = 217, we are actually in the range of case (5) with Ha 1/2 << Pe m << Ha 2, so that the left layer has split into a viscous layer with O(Ha - 1 / 2 ) thickness and a thinner mass-diffusion layer with O((HaPem ) - 1/3) = O ( y - J / 3 H a - 1 / 2 ) thickness. There is a more laterally uniform dopant concentration in the crystal, as shown in Fig. 9. The constant-concentration lines are only horizontal for 0 < x < 0.0015, and even in this initial region the maximum concentration at y = - 1 is only 7.9% higher than the C~ = 0.1029 at y = 1. After high concentrations begin returning to the left layer near v = 1, the crystal distribution becomes even more

0.93 0.92

0,15

y

0.0-

-0.5

0,12 0.5

jl

0.14

0.91 0.9 C, 0.89

-0,5 -

-1.0

0.88

)

-

L 0.00

I o.25

L

I o.5o

I

/ o.75

x

Fig. 8. Crystal composition C~(x y) for Bo = 4 T.

0.87 0.86

I

I

I

-0.5

O

0.5

Y

Fig. 10. Crystal Concentration C~(0.90 v) versus y for B0 = 0.5 T.

N. Ma, J.S. Walker/Journal of Crystal Growth 172 (1997) 124-135

laterally uniform, as shown by the lateral variation of C~ at x = 0.9 in Fig. 10. Here, the maximum at y=-I and Cs at y = - 0 . 5 are 7.0% and 1.6% larger than the C~ =0.8625 at y = 1. The axial variation in Fig, 9 is almost a perfect fit to the well-mixed formula, even though again the melt is far from well-mixed with extremely non-uniform concentrations in the melt at all times.

5. Conclusions Without turbulent mixing, convective mass transport dominates over diffusion, except in thin boundary layers with large concentration gradients. Convection leads to extremely non-uniform dopant distributions because the melt velocity is non-uniform. In our model problem for y < 0 , the fluid near y = - 1 quickly carries excess dopant across the bottom, leading to large values of C near y = - 1, while the fluid near y = - w moves very slowly, so that C here remains equal to one over most of the axial length of the melt. For y > 0 and for B 0 = 4 T, the concentration leaving the right layer at any time is relatively uniform, but it is increasing with time. Since the fluid near y = 1 moves much faster than that near y = 0, the C in the core for y > 0 and that entering the left layer is very non-uniform. With Pem >> 1, the concentration is constant along each streamline in the core for the instantaneous steadystate solution. The actual dopant distribution is always far away from the instantaneous steady state because particles moving along any streamline in the core exited from the left or right diffusion layers at different times and thus carry different concentrations. Therefore the dopant transport problem is intrinsically unsteady for magnetic fields which are strong enough to eliminate turbulence and laminar unsteadiness, but which are not strong enough to achieve diffusion-controlled mass transport. The pre-

135

sent asymptotic method can accurately predict the dopant distribution in an entire crystal with very modest computational resources.

Acknowledgements This research was supported by the National Aeronautics and Space Administration under Cooperative Research Agreement NCC8-90 and by the National Science Foundation under Grant CTS 9419484. The calculations were performed on the SGI Power Challenge supercomputer at the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign under Grant CTS 96-0024N.

References [1] H.A. Chedzey and D.T.J. Hurle, Nature 210 (1966) 933. [2] H.P. Utech and M.C. Flemings, J. Appl. Phys. 7 (1966) 2021. [3] D.H. Kim, P.M. Adornato and R.A. Brown, J. Crystal Growth 89 (1988) 339. [4] R.W. Series and D.T.J. Hurle, J. Crystal Growth 113 (1991) 305. [5] K.J. Lee, W.E. Langlois and K.M. Kim, Phys. Chem. Hydrodyn. 5 (1984) 135. [6] J.P. Garandet, J. Crystal Growth 125 (1992) 112. [7] T. Alboussi~re, Doctoral Thesis, Institut National Polytechnique de Grenoble (March 1992). [8] L.N. Hjellming and J.S. Walker, J. Fluid Mech. 182 (1987) 335. [9] J.S. Walker, G.S.S. Ludford and J.C.R. Hunt, J. Fluid Mech. 46 (1971) 657. [10] J.S. Walker, J. Mrcanique 20 (1981) 79. [11] T. Alboussi~re, J.P. Garandet and R. Moreau, J. Fluid Mech. 253 (1993) 545. [12] J.P. Garandet, J. Crystal Growth 114 (1991) 593. [13] N. Ma and J.S. Walker, Phys. Fluids 8 (1996) 944. [14] A. Ting, T.Q. Hua, J.S. Walker and B.F. Picologlou, Int. J. Eng. Sci. 31 (1993) 357.