Dopant segregation during liquid-encapsulated Czochralski crystal growth in a steady axial magnetic field

Journal of Crystal Growth 242 (2002) 471–485

Dopant segregation during liquid-encapsulated Czochralski crystal growth in a steady axial magnetic field Joseph L. Mortona, Nancy Mab,*, David F. Blissc, George G. Bryantc a

Department of Mechanical and Aerospace Engineering and Engineering Mechanics, University of Missouri, 1870 Miner Circle, Rolla, MO 65409, USA b Department of Mechanical and Aerospace Engineering, North Carolina State University, Campus Box 7910, Raleigh, NC 27695, USA c US Air Force Research Laboratory, AFRL/SNHC, 80 Scott Road, Hanscom AFB, MA 01731, USA Received 2 March 2002; accepted 1 May 2002 Communicated by D.T.J. Hurle

Abstract During the magnetically stabilized liquid-encapsulated Czochralski (MLEC) process, a single compound semiconductor crystal such as indium-phosphide (InP) or gallium-antimonide (GaSb) is grown by the solidification of an initially molten semiconductor (melt) contained in a crucible. The motion of the electrically conducting molten semiconductor can be controlled with an externally applied magnetic field. This paper presents a model for the unsteady transport of a dopant during the MLEC process with a steady axial magnetic field. The convective species transport during growth is driven by the melt motion, which produces segregation, i.e. non-uniformities in the dopant concentration, in both the melt and the crystal. This convective transport is significant even for a magnetic field strength of 2 T. Except for the last-solidified part of the crystal, the crystal’s axial dopant homogeneity, i.e. uniformity in the dopant concentration, improves as the magnetic field strength is decreased. Dopant distributions in the crystal and in the melt at several different stages during growth are presented for several magnetic field strengths. r 2002 Elsevier Science B.V. All rights reserved. PACS: 44.25; 81.10; 61.72.V; 47.65; 81.10.F; 64.75; 61.72; 72.20 Keywords: A1. Magnetic fields; A1. Mass transfer; A1. Segregation; A2. Growth from melt; A2. Liquid-encapsulated Czochralski method; A2. Magnetic field assisted Czochralski method; B2. Semiconducting III–V compounds

1. Introduction Since molten semiconductors are excellent electrical conductors, a magnetic field can be used to *Corresponding author. Tel.: +1-919-515-5231; fax: +1919-515-7968. E-mail address: nancy [email protected] (N. Ma).

control the melt motion in order to control the dopant distribution in the crystal, which depends on the convective and diffusive transport of the dopant in the melt. Convective transport in the melt may lead to (i) small-scale spatial oscillations of the crystal’s dopant composition, which are called striations or microsegregation, and/or (ii) large-scale variations of the crystal’s dopant

0022-0248/02/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 2 ) 0 1 4 2 5 - 2

472

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

composition, which are called radial and axial macrosegregation. Hurle and Series [1] and Walker [2] reviewed the literature on the use of magnetic fields during the bulk growth of semiconductor crystals. Walker and Ma [3] reviewed convective mass transport during semiconductor crystal growth with magnetic fields. During the bulk growth of semiconductor crystals without a magnetic field, the buoyant convection in the melt is often periodic, or even turbulent in very large systems. Unsteady melt motions lead to fluctuations in the heat transfer across the growth interface from the melt to the crystal. Since the local rate of crystallization depends on the balance between the local heat fluxes in the melt and crystal, fluctuations in the heat flux from the melt create fluctuations in the local growth rate. In extreme circumstances, periods of growth alternate with periods of remelting. Fluctuations in the local growth rate cause two major problems. First, these fluctuations are a major cause of dislocations in the crystal [4]. Another major cause of dislocations is the thermal stresses in the crystal. Second, the fluctuations in the local growth rate create striations in the crystal. Most dopants are either rejected into the melt during solidification or preferentially absorbed into the crystal, i.e. ks o1 or ks > 1; respectively, where the segregation coefficient ks is the ratio of the local dopant concentration in the crystal to that in the melt at any point along the crystal–melt interface. If ks o1; the rejected dopant accumulates in a species-diffusion boundary layer in the melt adjacent to the interface. The dopant distribution in this layer involves a balance between the rejection rate and the rate of diffusion through the melt. During a fluctuation in local growth rate, an increase in growth rate causes the local concentration to rise because dopant cannot diffuse away fast enough, so that a region with a high dopant concentration solidifies. When the growth rate decreases for the second half of the fluctuation, dopant has ample time to diffuse away, so the concentration that solidifies is low. Fluctuations in the melt velocity also convect dopant in and out of the species-diffusion boundary layer or radially within this layer, leading to more severe striations. It is not unusual for the

maximum dopant concentration to be twice the minimum in striations, and a typical distance between adjacent maximums is 1 mm. A magnetic field produced by a solenoid placed around the crystal-growth furnace can be used to stabilize the melt in order to eliminate all periodicity in the melt motion and to eliminate striations produced by unsteady melt motions. Without a magnetic field, transitions from steady, axisymmetric melt motions to the periodic or even turbulent non-axisymmetric melt motions, which produce striations, depend upon the ratio of the driving buoyancy force to viscous dissipation, as reflected by the Rayleigh or Grashof number. With a magnetic field, both viscous and Joulean dissipations oppose the driving buoyancy force and stabilize the flow through electromagnetic (EM) damping. The characteristic ratio of the Joulean dissipation to the viscous dissipation with an externally applied magnetic field is Ha2 ; where Ha ¼ B0 Rðs=mÞ1=2 is the Hartmann number, B0 is the magnetic flux density of the magnetic field, R is the characteristic dimension of the melt, s is the electrical conductivity of the melt and m is the dynamic viscosity of the melt. Since the ratio of s to m is enormous for all virtually molten semiconductors, even a weak magnetic field gives Ha ¼ 500 or 1000 for a typical crystal growth system, so that Lorentz force can easily stabilize the melt motion, in order to (i) prevent twinning, in which a second crystal grows with a different crystallographic orientation from that of the original crystal, (ii) dramatically reduce dislocation densities, (iii) eliminate striations in the crystal, and (iv) minimize macrosegregation in the crystal. Bliss et al. [5,6] were the first to produce 8-cm diameter twin-free indium-phosphide crystals by using magnetic stabilization. Hoffmann et al. [7] grew a nearly homogeneous crystal by using a strong magnetic field during the early stages of growth and reducing the field strength as crystal growth progressed. During directional solidification or Bridgman growth of semiconductor crystals, an extremely strong EM damping of the melt motion may produce a crystal with an axially and radially uniform dopant composition [8]. In order to achieve such growth conditions, the species

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

transport Pe! clet number, Pem ¼ Ub R=D; must be small, where Ub is the characteristic velocity for the magnetically damped melt motion, and is inversely proportional to the square of the magnetic flux density, while D is the diffusion coefficient for the dopant in the molten semiconductor. However, since typical values of D are 1–2  108 m2/s, B0 must be extremely large for diffusion-controlled species transport on Earth [9]. It is more practical to use a moderately strong magnetic field to tailor the melt motion in order to achieve both radial and axial dopant uniformity in the crystal during terrestrial crystal growth. 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, 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 that accurately predict the dopant distribution in an entire crystal for any combination of process variables are needed to facilitate process optimization. In the present paper, we treat the species transport of iron (Fe) in an indium-phosphide (InP) melt during the magnetically stabilized liquid-encapsulated Czochralski (MLEC) process with a uniform, steady, axial magnetic field. For a typical crystal growth process, resolution of thin species-diffusion boundary layers having O(Pe1 m ) thicknesses is often very challenging. Typically, Pem is very large for weak magnetic fields because the value of D is so small. Several grid points must be concentrated inside each layer in order to give accurate results, because these boundary layers play critical roles in the transport. As the magnetic flux density is increased, the value of Pem decreases and the species-diffusion boundary layers become thicker but the value of the Hartmann number increases. Since viscous boundary layers with O(Ha1 ) thicknesses exist, there are thin speciesdiffusion or viscous boundary layers for every value of B0 : Therefore, the simultaneous numerical solution of the full Navier–Stokes, internal energy and species transport equations must always have

473

a very fine spatial grid and a very small time step. Kaiser and Benz [10] stated that the hardware requirements and needed supercomputing time for the accurate resolution of concentration and velocity gradients are impractical. The present study eliminates the need for impractical computing resources by using an asymptotic approach to treat the unsteady dopant transport for the entire period needed to grow a crystal during the solidification of an InP crystal with a magnetic field. In a pair of studies, Ma and Walker [11,12] presented an asymptotic treatment which provided the first predictions of dopant composition in the entire crystal in an idealized configuration with a magnetic field. Recently, Ma and Walker [13] applied their method to a realistic crystal growth process, and provided the first predictions of dopant composition in the entire crystal grown by the vertical Bridgman process with a magnetic field. In the present paper, we apply this asymptotic approach to the MLEC process, and provide accurate predictions for the entire period needed to grow a crystal with moderate computational resources.

2. Temperature This paper treats the unsteady, axisymmetric species transport of iron in an indium-phosphide melt during the liquid-encapsulated Czochralski process with an externally applied, uniform, steady, axial magnetic field B0 z# : Here B0 is the magnetic flux density, while r#; h# and z# are the unit vectors for the cylindrical coordinate system. Experiments have demonstrated that a cusped axisymmetric magnetic field does not produce a better crystal [14]. Our dimensionless problem is sketched in Fig. 1. The coordinates and lengths are normalized by the crucible’s inner radius R; so that g is the dimensionless crystal radius, and bðtÞ is the dimensionless depth of the melt. This study uses the bulk approximation which assumes that the crystal–melt and encapsulant–melt interfaces lie in the same horizontal plane at z ¼ bðtÞ ¼ b0  og2 t; where t is normalized with R=Ub : Here, b0 is the initial dimensionless melt depth and the dimensionless crystal-growth velocity o ¼ Ug =Ub is the

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

474 quartz crucible

crystal

γ

B2O 3

B2O 3

Bo z z=b(t)

melt

z r

graphite susceptor

Fig. 1. MLEC crystal growth with a uniform, steady, axial magnetic field B0 z# and with coordinates normalized by the crucible’s inner radius.

sum of the velocity at which the crystal is moved upward and the velocity at which the crystal–melt interface moves downward. Before solidification begins, phosphorus gas is bubbled at high pressure through the indium melt, and the indium and phosphorus fuse to form the compound InP. A layer of boron oxide (B2O3) encapsulates the melt to prevent escape of the volatile component (P). A single crystal seed is lowered through the B2O3 which initiates solidification. Once the crystal has grown to the desired diameter, the crystal is pulled vertically upward at a rate which maintains this diameter. The melt is housed by a quartz crucible which is structurally supported by a graphite susceptor as shown in Fig. 1. For a typical process, ðDTÞc ¼ 20 K and R ¼ 4:7 cm [6]. In the energy equation, the characteristic ratio of the convective to conductive heat transfer is the thermal Pe! clet number, Pet ¼ rcp Ub R=k; where r is the melt’s density at the solidification temperature Ts ; cp is the melt’s specific heat, and k is the melt’s thermal conductivity. For a sufficiently strong magnetic field, Pet is small, so that the convective heat transfer terms are negligible. In a recent study, Ma and Walker [15] investigated the role of inertia and convective heat transfer on the buoyant convection during MLEC growth, and determined the errors associated with the neglect of convective heat transfer for field strengths between 0.1 and 0.9 T. For MLEC growth of indium phosphide, they found that convective heat

transfer significantly affect the buoyant convection when B0 ¼ 0:1 T for which Pet ¼ 277:3: As the magnetic field strength is increased from this value, the thermal Pe! clet number decreases as B2 0 ; so that the ratio of the convective heat transfer to conductive heat transfer decreases. Ma and Walker [15] found that the error due to neglect of convective heat transfer is o4% for B0 X0:43 T for which Pet p15:0: The present paper illustrates a method which can be used to predict the effect of buoyant convection on the dopant transport for growth in a 2- or 1-T magnetic field where convective heat transfer is totally negligible. Inclusion of convective heat transfer for weak magnetic field strengths is a straightforward extension of the present model, and will be investigated in a future study. In general, Ug oUb ; so that the latent heat released by the cooling melt is negligible compared to the conductive heat transfer [11]. The temperature is governed by r2 T ¼ 0;

ð1Þ

where T is the deviation of the dimensional temperature from the solidification temperature Ts normalized by the characteristic temperature difference in the melt ðDTÞc : We use the boundary conditions qT ¼ qb ðr; tÞ at z ¼ 1; qz

ð2aÞ

qT ¼1 qr

at r ¼ 1;

ð2bÞ

T ¼0

at z ¼ þ1

2 qT ¼ k0 bðtÞ qz

for 0prpg;

at z ¼ þ1

for gprp1;

ð2cÞ ð2dÞ

where qb ðr; tÞ is the dimensionless heat flux into the melt along the bottom crucible wall, and k0 is the dimensionless heat flux lost due to conduction and radiation through the semitransparent boron oxide. qb ðr; tÞ and k0 are normalized with ðDTÞc k=R: Here, z ¼ 1 þ 2z=bðtÞ is a rescaled axial coordinate so that 1pzp þ 1 for all time. Since the temperature variation along the melt–encapsulant interface is tiny compared

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

qb ðr; tÞ ¼ h½bðtÞ r2 ;

1 0.03 7

0.09 5

0

ζ

to the temperature variations in the furnace [16], we assume that the heat flux across the melt– encapsulant interface is uniform. An estimate of the thermal conduction through the encapsulant based on measured temperatures [16] and an estimate of the radiative losses through the semitransparent encapsulant [17] give k0 ¼ 1:2: During the MLEC process, radio-frequency (RF) induction heating is used, and has the advantage that the heating distribution can be tailored by adjusting the number of loops and the spacing between the loops in the induction coil [18]. With RF heating, the distribution of heat flux into the horizontal and vertical susceptor walls is highly non-uniform, so that the distribution of the heat flux which has conducted through the susceptor, and conducted and radiated through the quartz crucible and into the bottom and side of the melt are also non-uniform. In the present study, we treat the thermal problem with uniform side heating and with parabolically varying bottom heating which varies over time so that 0pTðr; z; tÞp1 for all stages of growth. For MLEC growth, the bottom heating is estimated by a quadratic form given by

475

0.2 0.4 0.6 0.8

−1 0

1 r

Fig. 2. Isotherms for g ¼ 0:4 and b=0.3174.

3. Melt motion The melt velocity is normalized by the characteristic velocity for magnetically damped buoyant convection [19]:

ð3aÞ Ub ¼

where h½bðtÞ ¼ 0:890b þ 0:369b2  5:984b3 4

5

þ 7:334b  3:273b :

ð3bÞ

Future research will further investigate the effects of RF heating distribution on the buoyant convection and on the dopant transport with the objective to tailor the heating distribution in order to achieve a better crystal, e.g., one with a more uniform dopant distribution. A Chebyshev spectral collocation method is used to solve for the temperature governed by Eqs. (1) and (2) with Gauss–Lobatto collocation points in r and z: A regularization method is implemented to avoid a Gibb’s phenomena associated with the discontinuous boundary condition at r ¼ g and z ¼ þ1: Typical isotherms are presented in Fig. 2 for g ¼ 0:4 and b ¼ 0:3174; which correspond to 50% of growth.

rgbðDTÞc ; sB20

ð4Þ

where g ¼ 9:81 m/s2 and b is the thermal volumetric expansion coefficient. The crystal–melt interface moves at a constant dimensionless velocity o; and the dimensionless time to grow the entire crystal is b0 =og2 if the entire melt is solidified. We have assumed that there is no change in density during solidification. In the present study, we treat only two flows, namely the buoyant convection and the so-called melt-depletion flow. In a reference frame moving with the crucible, the fluid at the encapsulant–melt interface moves downward at the rate db=dt while the fluid at the crystal–melt interface moves vertically upward at the pull velocity [20], causing a melt motion which is referred to as the meltdepletion flow. This terminology arises because the crystal–melt interface acts as a porous boundary with suction.

476

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

Unfortunately, these two flows produce radially inward melt motion below the crystal–melt interface which provides a convection of dopant which produces severe radial segregation. During growth, the crystal and crucible are typically rotated for axisymmetry so we investigated the balance of the melt motion due to these rotations and buoyant convection for MLEC [21]. Because the magnitude of the centrifugal pumping varies as B3 while the magnitude of the buoyant 0 convection varies as B2 0 ; a rotation rate of 80 rpm is needed to produce sufficient radial outward flow for B0 ¼ 0:5 T [21]. This large rotation rate is not practical. With melt-depletion flow, an even larger rotation rate would be needed. Therefore, while the crystal and crucible rotations are needed for axisymmetry, other promising means are needed to improve radial dopant homogeneity. The electric current in the melt produces an induced magnetic field which is superimposed upon the applied magnetic field produced by the external magnet. The characteristic ratio of the induced to applied magnetic field strengths is the magnetic Reynolds number, Rm ¼ mp sUb R; where mp is the magnetic permeability of the melt. For all crystal growth processes, Rm 51 and the additional magnetic fields produced by the electric currents in the melt are negligible. In the Navier–Stokes equation, the characteristic ratio of the EM body force term to the inertial terms is the interaction parameter, N ¼ sB20 R= rUb ; which varies as B40 : For a sufficiently strong magnetic field, N is large, so that the inertial terms are negligible. In a recent study, Ma and Walker [15] investigated the role of inertia during MLEC growth, and determined the errors associated with the neglect of inertial effects for field strengths between 0.1 and 0.9 T. They found that inertia significantly affects the buoyant convection when B0 ¼ 0:1 T for which N ¼ 1:037: As the magnetic field strength is increased from this value, the interaction parameter increases as B40 ; so that the ratio of the inertial force to the EM body force decreases. Ma and Walker [15] found that the error due to neglect of inertial effects is only 2.7% for B0 ¼ 0:2 T where N ¼ 16:59; and is totally negligible for B0 X0:5 T.

For an inertialess flow, the equations governing the melt motion are: rp þ fb z# þ j  z# þ Ha2 r2 v ¼ 0;

ð5aÞ

r v ¼ 0;

ð5bÞ

r j ¼ 0;

ð5cÞ

j ¼ rf þ v  z# ;

ð5dÞ

where vðr; z; tÞ ¼ vr r# þ vz z# is the dimensionless velocity of the melt normalized by Ub ; p is the deviation of the dimensional pressure from the hydrostatic pressure normalized by sB20 Ub R; j is the electric current density normalized by sUb B0 ; and f is the electric potential function normalized by Ub B0 R: In the Navier–Stokes equation, f b is a body force which is normalized by sUb B20 : fb ¼ 0 for melt-depletion flow, and fb ¼ T for buoyant convection with the Boussinesq approximation, where T is given by a solution to Eqs. (1) and (2). In an asymptotic solution for Hab1; the melt is divided into an inner inviscid core region below the crystal–melt interface, an outer inviscid core region below the encapsulant–melt interface, an O(Ha1=2 ) free shear layer at r ¼ g beneath the crystal’s edge, Hartmann layers with an O(Ha1 ) thickness adjacent to the bottom crucible wall at z ¼ 1 and adjacent to the crystal–melt and encapsulant–melt interfaces at z ¼ þ1; and a parallel layer with an O(Ha1=2 ) thickness adjacent to the crucible’s vertical wall. For most field strengths of interest, the encapsulant is so viscous that vr ¼ 0 at the encapsulant–melt interface, so that we treat the encapsulant–melt interface as a no-slip boundary [22]. A strong EM damping of the melt produces a motion in the melt that is comparable to that of the encapsulant so that the flows are coupled for very strong fields [22]. The Hartmann layers have a simple, local, exponential structure, which matches any radial core, free shear layer or parallel layer velocities at z ¼ 71; which satisfies the no-slip conditions at the crucible’s surface, the crystal–melt interface and the encapsulant–melt interface, and which indicates that vz in the cores, free shear layer or parallel layer is O(Ha1 ) at z ¼ 71: Since the core

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

carries an O(1) flow, the flow circuit must be completed through a high-velocity parallel layer with vz ¼ OðHa1=2 Þ and vr ¼ Oð1Þ inside this layer at r ¼ 1: For magnetic field strengths of interest, this parallel layer is not actually thin as assumed in the formal asymptotic expansion for Hab1: While a formal asymptotic analysis for Hab1 is not appropriate, the numerical solution of the inertialess Navier–Stokes equation with all viscous terms is not necessary. The Hartmann layers represent an extremely small fraction of the melt depth and have a simple exponential structure. There is no need to numerically duplicate this simple exponential structure. Therefore, we use a hybrid solution which does not assume that the parallel layer thickness is small. We discard the viscous terms Ha22 v=z2 in the Navier–Stokes equation, and we relax the no-slip conditions at z ¼ 71 because they are satisfied by the Hartmann layers which are not part of the hybrid solution. For our axisymmetric flow, there is no azimuthal electric field, and there is a radial EM body force due to the azimuthal electric current and the axial magnetic field. The streamfunction c is governed by     q 1 qc 4 1 q2 c qfb 2 q 1 q r Ha ; ð6aÞ ¼  2 qr r qr qr r qr b r qz2 qr where vr ¼

2 1 qc ; b r qz

vz ¼ 

1 qc : r qr

ð6b; 6cÞ

The boundary conditions along the crystal–melt and encapsulant–melt interfaces are: c ¼ 12 oð1  g2 Þr2

at z ¼ þ1

for 0prpg; ð7aÞ

c ¼ 12 og2 ð1  r2 Þ

at z ¼ þ1

for gprp1: ð7bÞ

We use a Chebyshev spectral collocation method to solve Eq. (6a) with Eqs. (7a) and (7b) at z ¼ þ1; the no-slip and no-penetration conditions along the crucible’s vertical wall, and the no-penetration condition along the crucible’s bottom wall. We use Gauss–Lobatto collocation points in r and z; and a sufficient number of collocation points in order to resolve the velocity gradients near z ¼ þ1: Without buoyant convection, the melt-depletion flow alone is governed by Eq. (6a) with fb ¼ 0 and

477

boundary conditions (8a) and (8b). For a constant growth rate, the minimum value of the streamfunction is constant for all stages of growth. The melt-depletion flow is spread over almost the entire volume of the melt with axially downward flow for gprp1; radially inward flow everywhere, and axially upward flow for 0prpg: Without a magnetic field, this flow would be confined to a region which is much more local to the crystal– melt and encapsulant–melt interfaces. When the electrically conducting fluid flows radially across the vertical magnetic field, it generates an induced electric field which drives an azimuthal electric current which in turn flows across the magnetic field, creating an EM body force which opposes the radial velocity. Since there is no EM body force opposing flow along magnetic field lines, the flow along the magnetic field lines is elongated, so that there is melt-depletion flow over the entire melt. For the MLEC growth of indium-phosphide crystals, Ug ¼ 5:556 mm/s and all results use this value of the crystal-growth velocity. For g ¼ 0:4; b ¼ 0:3174 and Ha ¼ 2; 748; the minimum value of the melt-depletion streamfunction is cmin ¼ 0:00235623: Without the melt-depletion flow, the buoyant convection alone is governed by fb ¼ T with c ¼ 0 along z ¼ þ1 for 0prp1: Because the magnitude of the buoyant convection decreases with the melt depth, the maximum value of the streamfunction decreases as crystal growth progresses. The hot fluid rises near the periphery of the crucible, flows radially inward adjacent to the encapsulant–melt interface, and continues to flow radially inward near the crystal–melt interface where it is cooled. The fluid adjacent to the crystal–melt interface either solidifies into the crystal or flows radially downward, and then flows radially outward adjacent to the crucible’s bottom wall. The balance between the melt-depletion flow and the buoyant convection depends on the crystal-growth velocity, the magnetic field strength and the melt depth. Typically Ub > Ug ; so that during early stages of growth, the buoyant convection is much stronger than the melt-depletion flow, and the melt motion is dominated by the buoyant convection. For example, the streamlines for an initial melt depth b0 ¼ 0:6348 and g ¼ 0:4;

478

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

Ha ¼ 2; 748; and Ug ¼ 5:556 mm/s very nearly resemble those of the buoyant convection alone, and have a maximum streamfunction value cmax ¼ 0:0489246: The only noticeable difference in the streamlines is in a small region below the crystal– melt interface where the buoyant convection is small, and there is vertically upward flow due to the melt-depletion flow. This phenomena is

reflected in Fig. 3a at t ¼ 5:608 when 5% of the crystal has grown, for which cmax ¼ 0:0442069: For a given magnetic field strength, the magnitude of the buoyant convection velocity is proportional to bðtÞ and decreases as crystal growth progresses while the magnitude of the melt-depletion flow is proportional to Ug and is constant throughout growth. By the time 50% of the crystal has grown,

Fig. 3. Streamlines for the melt motion with Ha ¼ 2748; g ¼ 0:4 and Ug ¼ 5:556 mm/s: (a) cðr; z; 5:608Þ; (b) cðr; z; 56:08Þ; and (c) cðr; z; 100:94Þ:

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

the melt depth has decreased to b ¼ 0:3174; so that there is noticeable melt-depletion flow over a larger volume of the melt below the crystal–melt interface as reflected in Fig. 3b where the maximum value of the streamfunction is cmax ¼ 0:0128973: The melt flows radially outward along z ¼ 1; vertically upward along r ¼ 1; radially inward along z ¼ þ1 adjacent to the encapsulant– melt interface, and then either solidifies or flows vertically downward. At this time, the vertically upward melt-depletion flow has become stronger than the vertically downward velocity of the colder sinking fluid. For r > 0:5; the velocity of the rising fluid adjacent to r ¼ 1 has decreased due to the opposing melt-depletion flow. However, for 0:5prpr0 where r0 is the radial position corresponding to cmax ; the velocity of the sinking melt and the vertically downward melt-depletion flow are superimposed to produce a downward flow larger than that produced by either the buoyant convection or melt-depletion flow alone. By the time 80% of the crystal has grown, the maximum value of the streamfunction has decreased to cmax ¼ 0:00226413; and the motion due to the melt-depletion flow represents more than half the volume of the melt. When t ¼ 100:94; 90% of the crystal has grown and the streamlines are presented in Fig. 3c, where the maximum value of the streamfunction is cmax ¼ 0:000572207: The minimum value of the streamfunction remains at the melt-depletion flow’s value of cmin ¼ 0:00235623 for the entire process. Consequently, there is convective species transport for all stages of growth, unlike directional solidification or Bridgman crystal growth where there is no meltdepletion flow. For Bridgman crystal growth, the species transport becomes dominated by diffusive species transport near the end of growth because the buoyant convection becomes small [12].

4. Dopant transport Before solidification begins, the dopant concentration is uniform, and this initial value is used to normalize the concentration C; so that Cðr; z; t ¼ 0Þ ¼ 1: Once crystal growth starts, the dopant is rejected into the melt for ks o1; leaving a dopant-

479

rich region near the crystal–melt interface. The flow convects the melt with the high Fe concentration, thus referred to as ‘‘dopant-rich’’, away from the interface. The dimensionless equation governing this dilute species transport is qC 2 þ v rC ¼ Pe1 m r C; qt

ð8Þ

where v is dimensionless velocity for the net melt motion which is equal to the sum of the dimensionless velocities for the melt-depletion flow and for the buoyant convection. Since the streamfunctions for the buoyant convection and melt-depletion flow are positive and negative, respectively, there is some cancellation when they are summed. The boundary conditions for the crystal–melt and encapsulant–melt interfaces are 2 qC ¼ Peg ð1  kS ÞC b qz

at z ¼ þ1 for 0prpg; ð9Þ

where Peg ¼ Ug R=D ¼ oPem is the growth Pe! clet number. The boundary conditions along the encapsulant–melt interface and along the bottom and vertical crucible walls are n# rC ¼ 0; where n# is the outward unit normal vector. Solutions for the dopant concentrations in the Hartmann layers show that the changes in the boundary conditions along z ¼ 71 due to the Hartmann layers are smaller than the O(Ha1 ) terms which are neglected in the melt motion solutions. Therefore, the Hartmann layers can be ignored in the species transport problem. We use a Chebyshev spectral collocation method for the convective and diffusive terms in Eq. (8) with Gauss–Lobatto collocation points in r and z: We use a sufficient number of collocation points in order to resolve the concentration gradients near z ¼ 71 and the velocity gradients near r ¼ 1: For the time derivative in Eq. (8), we use a secondorder implicit time integration scheme from t ¼ 0 to a t which is slightly less than b0 =og2 : We use a large enough number of times steps such that the results do not change by increasing the number of time steps. At the beginning of crystal growth, the melt concentration, normalized with the initial uniform

480

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

concentration, is Cðr; z; t ¼ 0Þ ¼ 1: Thus, the amount of dopant initially in the melt is obtained by integrating across the melt’s volume giving a total dopant concentration equal to pb0 : We verify that the sum of the dopant concentration in the melt and in the crystal is equal to pb0 at each time step. Assuming that there is no diffusion of dopant in the solid crystal and that the densities of the solid and liquid are the same, the dopant distribution in the crystal Cs ðr; XÞ; normalized by the initial uniform dopant concentration in the melt, is given by   b0  g 2 X CS ðr; XÞ ¼ kS C r; z ¼ þ1; t ¼ ; ð10Þ og2 where ks ¼ 0:001 for indium-phosphide doped with iron [23]. Here, X ¼ b0 =g2 corresponds to the first-grown part of the crystal at t ¼ 0 while X ¼ 0 corresponds to the last-grown part of the crystal. The first-grown part of the crystal at t ¼ 0 solidifies with Cs ðr; b0 =g2 Þ ¼ ks : In the following sections, we present results for growth of Fe-doped InP by the MLEC process for which g ¼ 0:4; b0 ¼ 0:6348; ks ¼ 0:001 and Ug ¼ 5:556 mm/s, for which Peg ¼ 13:0556: We present results for several cases in order to investigate the effect of the magnetic field strength B0 : The parameters as a function of B0 are Ub ¼ 0:000628242B2 0 m/s, Ha ¼ 1374:1B0 ; and Pem ¼ 1476:4B2 while the dimensionless crystal0 growth velocity is o ¼ 0:008843B20 ; and the dimensionless time to grow the crystal is b0 =og2 ¼ 448:659B2 0 ; where B0 is in Tesla.

5. Results for B0 ¼ 2 T We present results for B0 ¼ 2 T, for which o ¼ 0:035372; Ha ¼ 2; 748 and Pem ¼ 369:1: The characteristic velocity for the magnetically damped buoyant convection is Ub ¼ 0:000157 m/s, and the dimensionless time to grow the crystal is b0 =og2 ¼ 112:16: The constant–concentration curves in the melt and in the crystal are presented in Figs. 4 and 5, respectively. When 5% of the crystal has grown, most of the melt below the encapsulant–melt interface remains at the initial uniform concentration C ¼ 1 and the

maximum value of the concentration is 6.828, as shown in Fig. 4a. The dopant-rich melt has already diffused axially downward past a small region below the crystal–melt interface where the melt-depletion flow dominates. The axially downward diffusion is opposed by an axially upward velocity. In Fig. 4a, the dopant which has escaped from this region below the crystal–melt interface has been convected axially downward by the buoyant convection as reflected in the C ¼ 1:54 and 1:16 contours and has been further convected radially outward as reflected in the C ¼ 1:03; 1:01 and 1:001 contours. Near the edge of the crystal, the dopant-rich melt diffuses radially outward and downward. This diffusion is opposed by a relatively weak axial velocity which is roughly one-fourth the velocity at the centerline, so that dopant is easily transported out of this region. Furthermore, once this dopant has diffused a short distance downward, it is easily convected downward with a large downward velocity. The radially inward flow adjacent to the interfaces contributes to the accumulation of dopant below the crystal– melt interface. When 9% of the crystal has grown, the dopant has convected or diffused so that C > 1 everywhere. When 50% of the melt is remaining at t ¼ 56:08; the minimum and maximum values of the concentration are 1.559 and 21.94, respectively, as shown in Fig. 4b. The buoyant convection carries the dopant-rich melt into all regions of the melt. The upward buoyant convection adjacent to the crucible’s vertical wall convects the dopantrich melt towards the encapsulant–melt interface as reflected in the C ¼ 1:567; 1:578; 1:595 and 1:62 contours in Fig. 4b. When 90% of the crystal has grown, the minimum and maximum values of the concentration are 2.322 and 286.9, respectively, as shown in Fig. 4c. The maximum concentration is very large because of the extremely small value of the segregation or partition coefficient. The axially upward flow below the crystal–melt interface provides dopant convection that opposes diffusion away from the interface. This has caused the constant–concentration curves to appear nearly vertical for most of the melt except in a thin region below the crystal–melt interface.

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

481

Fig. 4. Contours of concentration in the melt Cðr; z; tÞ with Ha ¼ 2748; g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/s and ks ¼ 0:001: (a) Cðr; z; 5:608Þ; (b) Cðr; z; 56:08Þ; and (c) Cðr; z; 100:94Þ:

For larger values of X which correspond to the earlier grown sections of the crystal in Fig. 5, the constant–concentration curves have a large downward slope. This can be attributed to the large inward radial velocity opposing the radially outward diffusion in the top of the melt during the early stages of growth. Near the periphery of the crystal, the

concentration remains near ks for a long decrease in X as reflected by the Cs ¼ 0:0016 contour. The difference in the crystal’s composition at its centerline and periphery, GðXÞ; is summarized in Table 1 for several axial positions in the crystal. The radial segregation reflected by this difference GðXÞ becomes higher as crystal growth progresses.

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

482

1 2.68

0. 0 0 1 6 3

1.01 7

1.6

0 .0029

1.015

1.1 0. 0 0 6 4

1.01 3

Ξ

0

ζ

2

1.01 1 0. 0 1 6 3

1.00 9 0. 0 3 9 4

1

0. 12 7

−1 0

0

0.2 r

1

r

0.4

Fig. 5. Contours of concentration in the crystal Cs ðr; XÞ with Ha ¼ 2748; g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/s, and ks ¼ 0:001:

Fig. 6. Contours of concentration in the melt Cðr; z; t ¼ 22:433Þ with Ha ¼ 1374; g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/ s, and ks ¼ 0:001:

0. 00 12

Table 1 Radial variation of the dopant composition in the crystal GðXÞ ¼ Cs ð0:4; XÞ  Cs ð0; XÞ B0 ¼ 1 T

0.005221 0.01162 0.01917 0.07465 0.2804

0.003995 0.005656 0.01017 0.04963 0.2186

0. 00 2 2

Ξ

Gð3:7691Þ Gð2:9756Þ Gð1:98375Þ Gð0:7935Þ Gð0:39675Þ

B0 ¼ 2 T

0. 0 0 15 3

0 .0036

2

0. 0 0 67 0. 0 2

1

6. Results for B0 ¼ 1 T We decrease the field strength to 1 T, for which Ha ¼ 1; 374 and Pem ¼ 1476:4; so that the characteristic velocity has increased to Ub ¼ 0:000628 m/s, and the dimensionless time to grow the crystal is b0 =og2 ¼ 448:66: The constant– concentration curves in the melt and in the crystal are presented in Figs. 6 and 7, respectively. With B0 ¼ 1 T, the buoyant convection is much stronger than the melt-depletion flow for most of growth so that most of the melt flows in a counterclockwise circulation. The radially inward and axially upward melt-depletion flow occurs in a

0. 06 0. 27 5

0

0.2

0.4

r Fig. 7. Contours of concentration in the crystal Cs ðr; XÞ with Ha ¼ 1374; g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/s, and ks ¼ 0:001:

very small region below the crystal–melt interface. When 5% of the crystal has grown, the maximum value of the streamfunction is cmax ¼ 0:0435889:

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

When 50% of the crystal has grown, the maximum value of the streamfunction is cmax ¼ 0:01295 and the buoyant convection continues to dominate the flow except in an extremely small region below the interfaces. Even when 70% of the crystal has grown, the flow below the crystal–melt interface is axially upward only for z > 0:5: When 90% of the crystal has grown, the maximum value of the streamfunction is cmax ¼ 0:0007759 and the buoyant convection’s counterclockwise circulation extends in the range 0:5prp1: At the equivalent stage during growth for B0 ¼ 2 T, this circulation takes much less than half the melt’s volume. The minimum value of the streamfunction is cmin ¼ 0:0005891 for all stages of growth. With a 1 T field strength, the buoyant convection provides a much stronger convection of dopant so that the dopant has already convected or diffused so that C > 1 everywhere in the melt when only 4% of the crystal has grown. The constant–concentration curves at 5% of growth when t ¼ 22:433 are presented in Fig. 6. In Fig. 6, the minimum and maximum values of the concentration are 1.007 and 5.398, respectively, and the buoyant convection has carried the dopant to the entire volume of the melt so that there are closed constant–concentration curves. The radially inward flow for z > 0 opposes outward radial diffusion adjacent to the crystal–melt interface, while the outward radial flow for zo0 is consistent with the direction of the outward radial diffusion. This is reflected in the shape of the C ¼ 1:1 contour in Fig. 6, which curves radially inward as z decreases from z ¼ þ1; becomes more vertical near z ¼ 0; and then curves radially outward near the bottom of the crucible. The buoyant convection has convected the dopant so that many of the contours have closed on themselves below the encapsulant–melt interface as shown in Fig. 6. In fact, the convection is so strong that we see a local hump in these contours due to the large vertically upward velocity in the parallel layer adjacent to the vertical crucible wall. Even though the melt depth has decreased by 50% at t ¼ 224:33; the buoyant convection continues to dominate. At this stage, the minimum and maximum concentrations in the melt are 1.805 and 12.51, respectively. With the increased mixing,

483

the dopant is not being trapped below the crystal– melt interface as much as for growth with the stronger magnetic fields. In fact, this maximum value is nearly half the value for B0 ¼ 2 T at the same stage of growth. As crystal growth progresses, the decreasing strength of the buoyant convection allows a stronger transport of dopant due to the melt-depletion flow. When 90% of the crystal has grown, the minimum and maximum concentrations in the melt are 3.6250 and 225.8, respectively, and the constant–concentration curves resemble those of Fig. 4b. This difference between the maximum and minimum concentrations is less than the difference for the stronger field strengths due to the increased mixing. The constant–concentration curves in the crystal are presented in Fig. 7. Both the radial and axial variation of the concentration in the crystal are significantly different than the variations for B0 ¼ 2 T. While the dopant concentration near the centerline increases as X decreases, the concentration increases at a smaller rate than in the previous case because the increased convection carries dopant away from z ¼ þ1 at the centerline in the melt. Near the edge of the crystal, the dopant concentration remains near ks for a longer decrease in X: Although the concentration levels near r ¼ 0 do not rise as dramatically as in the prior case, the radial segregation for the earlier grown sections of the crystal remains large, as shown in Table 1. However, there is less radial segregation for the present case with B0 ¼ 1 T than with B0 ¼ 2 T.

7. Comparison with classical limits The classical limits where the dopant transport is dominated by either diffusive mass transport or convective mass transport are usually defined in terms of directional solidification or Bridgman crystal growth, where there is no liquid encapsulant [24]. For directional solidification with a planar interface, there is no radial variation of the crystal’s composition for diffusion-controlled growth. This limit has a different meaning for the liquid-encapsulated Czochralski process, where the crystal–melt interface and encapsulant–melt

484

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

interface lie in the same horizontal plane. Because of the presence of the melt-depletion flow, even in the absence of forced or natural convection, convective dopant transport creates both an axial and a radial variation of the dopant composition everywhere in the crystal. The axial variation of the crystal’s radially averaged composition for the Czochralski equivalent of the diffusion-controlled limit is presented in curve a in Fig. 8 for Ha ¼ 5; 496; g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/s and ks ¼ 0:001: There is a significant axial macrosegregation, which is very different from directional solidification which would produce an axially uniform crystal except in the first-grown and last-grown parts of the crystal [24]. During crystal growth without a magnetic field, the mixing may produce a dopant composition in the crystal that is close to the classical well-mixed limit. If the rejected dopant is instantly uniformly distributed over the entire volume of the melt at each stage during growth, then the crystal’s dopant composition is radially uniform and the axial variation is   b0 ð1ks Þ Cs;well-mixed ðXÞ ¼ ks 2 : ð11Þ g X

The axial variation of the crystal’s composition for the well-mixed limit is presented in curve e in Fig. 8 for g ¼ 0:4; b0 ¼ 0:6348 and ks ¼ 0:001: As the magnetic field strength is increased from zero, the buoyant convection is reduced and approaches the B2 variation in Ub ; so that the 0 dopant transport transitions from the well-mixed limit for zero or weak magnetic fields to the diffusion-controlled limit for extremely strong magnetic fields. The axial variation of the radially averaged crystal composition for various magnetic field strengths are presented in Fig. 8. For each case, the crystal initially solidifies with Cs ¼ ks ¼ 0:001 at X ¼ 3:9675: Curve a in Fig. 8 corresponds to diffusion-controlled growth for B0 ¼ 4 T. The maximum concentration in the melt at each stage during growth exists at the centerline at z ¼ þ1; so that the crystal solidifies with an elevated C at the centerline for every value of X: The maximum concentration in the melt at each stage during growth decreases for decreasing B0 ; so that the crystal solidifies with consistently lower concentration for a given X for decreasing B0 except near the end of growth. This phenomena is shown in curves b, c and d in Fig. 8 which correspond to growth with buoyant convection for B0 =4, 2 and 1 T, respectively. For ks ¼ 0:001; the crystal absorbs very little dopant at each stage during growth so that the dopant concentrations in the melt are extremely large near the end of growth. At the end of growth, all of the remaining dopant-rich melt must solidify, so that the last-grown sections of the crystal solidify with much higher concentrations, as reflected in the steepness of the curves in Fig. 8 for small values of X:

8. Conclusions

Fig. 8. Axial variation of the radially averaged crystal concentration Cs;avg ðr; XÞ with g ¼ 0:4; b0 ¼ 0:6348; Ug ¼ 5:556 mm/s and ks ¼ 0:001: Cases: (a) diffusion-controlled limit with Ha ¼ 5496; (b) Ha ¼ 5496; (c) Ha ¼ 2748; (d) Ha ¼ 1374; and (e) well-mixed limit.

We have developed a method which we can use to accurately predict dopant transport during liquid-encapsulated crystal growth in a steady magnetic field. We found that crystal growth in a 2 T magnetic field causes a deviation of the crystal composition from radial uniformity. As the magnetic field strength is decreased, the radial and axial homogeneity of the crystal generally improves. For the crystal grown in 2 T with or

J.L. Morton et al. / Journal of Crystal Growth 242 (2002) 471–485

without buoyant convection, there is significant segregation. As the field strength is decreased, there is less electromagnetic damping and the strength of the buoyant convection increases, so that some of the dopant that would otherwise accumulate below the crystal–melt interface is convected away by the buoyant convection. However, the extremely small value of ks creates extremely elevated dopant concentrations near the end of growth in the melt and in the last-grown sections of the crystal for all cases. Except near the ends of the crystal, there is less radial and axial segregation for the crystal grown in the weakest magnetic field considered.

Acknowledgements This research was supported by the National Aeronautics and Space Administration under Grant NAG8-1817, and by the US Air Force Office of Scientific Research. The calculations were performed on the IBM SP at the North Carolina Supercomputing Center in Research Triangle Park, NC. The authors are grateful to Professor J.S. Walker in the Department of Mechanical and Industrial Engineering at the University of Illinois at Urbana-Champaign for his thoughtful and constructive suggestions.

References [1] D.T.J. Hurle, R.W. Series, Handbook of Crystal Growth, Elsevier Science Publishers, Amsterdam, Vol. 2A, 1994, p. 261. [2] J.S. Walker, in: K.W. Benz (Ed.), The Role of Magnetic Fields in Crystal Growth, Progress in Crystal Growth and Characterization of Materials, Vol. 38, Elsevier, Amsterdam, 1999, p. 195.

485

[3] J.S. Walker, N. Ma, in: Chang-Lin Tien, Vishwanath Prasad, Frank Incropera (Eds.), Annual Review of Heat Transfer, Begell House, New York, Vol. 12, 2002. pp. 223–263. [4] E. Kuroda, H. Kozuka, Y. Takano, J. Crystal Growth 68 (1984) 613. [5] D.F. Bliss, R.M. Hilton, S. Bachowski, J.A. Adamski, J. Electron. Mater. 20 (1991) 967. [6] D.F. Bliss, R.M. Hilton, J.A. Adamski, J. Crystal Growth 128 (1993) 451. [7] D. Hoffmann, D. Mosel, G. Muller, in: D. Grossman, L. Ledebo (Eds.), Semi-insulating III–V Materials, IOP Publishing, Bristol, 1988, p. 429. [8] D.H. Kim, P.M. Adornato, R.A. Brown, J. Crystal Growth 89 (1988) 339. [9] N. Ma, J.S. Walker, J. Thermophys. Heat Transfer 11 (1997) 212. [10] Th. Kaiser, K.W. Benz, J. Crystal Growth 183 (1998) 564. [11] N. Ma, J.S. Walker, J. Crystal Growth 172 (1997) 124. [12] N. Ma, J.S. Walker, J. Crystal Growth 180 (1997) 401. [13] N. Ma, J.S. Walker, J. Heat Transfer 122 (2000) 159. [14] G.G. Bryant, D.F. Bliss, D. Leahy, R. Lancto, N. Ma, J.S. Walker, IEEE Proceedings of the International Conference on Indium Phosphide and Related Materials, 1997, pp. 416–419. [15] N. Ma, J.S. Walker, J. Thermophys. Heat Transfer 15 (2001) 50. [16] V. Prasad, D.F. Bliss, J.A. Adamski, J. Crystal Growth 142 (1994) 21. [17] F. Dupret, P. Nicod"eme, Y. Ryckmans, M.J. Crochet, Int. J. Heat Mass Transfer 33 (1990) 1849. [18] N. Ma, A.P. Anselmo, D.F. Bliss, J.S. Walker, ASME Proceedings of the 31st National Heat Transfer Conference, Vol. I, Thermal Transport in Materials Processing, HTD-Vol. 323, 1996, pp. 309–313. [19] L.N. Hjellming, J.S. Walker, J. Fluid Mech. 182 (1987) 335. [20] N. Ma, J.S. Walker, D.F. Bliss, G.G. Bryant, J. Fluids Eng. 120 (1998) 844. [21] J.L. Morton, N. Ma, D.F. Bliss, G.G. Bryant, J. Fluids Eng. 123 (2001) 893. [22] M.V. Farrell, N. Ma, J. Heat Transfer 124 (2002). [23] D.T.J. Hurle, Crystal Pulling from the Melt, Springer, New York, 1993. [24] M.C. Flemings, Solidification Processing, McGraw-Hill Inc., New York, 1974.