Effects of Rotation on Heat Flow, Segregation, and Zone Shape in a Small-scale Floating-zone Silicon Growth under Axial and Transversal Magnetic Fields

Effects of Rotation on Heat Flow, Segregation, and Zone Shape in a Small-scale

Floating-zone Silicon Growth under Axial and Transversal Magnetic Fields

C. W. Lan1_{, B. C. Yeh}

Abstract: The suppression of unstable Marangoni convection in floating-zone crystal growth by magnetic fields has enjoyed over recent years a widespread use as a reliable and useful strategy. A transversal direction of the field is particularly efficient, but asymmetric zone shapes and thus segregation are induced. Counter-rotation of the feed and of the crystal rods is a common way to improve dopant homogeneity. However, its effects under mag-netic fields are complex and have not yet been studied in detail. In the present analysis, three-dimensional (3D) simulations based on a finite-volume/multigrid method are used to illustrate the effects of rotation on the heat flow, dopant segregation, and the zone shape for a small-scale floating-zone silicon growth under both axial and transversal magnetic fields. The role of electrical con-ductivity of the crystal is also taken into account.

keyword: Rotation, Convection, Segregation, Mag-netic field, Lorentz force, Floating-zone method.

1 Introduction

The floating-zone (FZ) technique is a crucible-less and thus a contamination-free process for growing single crystals. However, the free surface of the molten zone often induces significant Marangoni flow leading to un-stable convection, growth striations, and distorted inter-faces, even in microgravity environment (Schweizer et al., 1999). Suppressing the unsteady Marangoni flow has been an important research topic in crystal growth over the years (Nakamura et al., 1998; Amberg and Shiomi, 2005), and static magnetic fields seem to be particularly useful in this regard. Striation-free crystals can be grown under a strong static magnetic field (see Schweizer et al., 1999; Dold et al., 1998).

In practice, both axial (Schweizer et al., 1999; Dold et

1_{Department of Chemical Engineering National Taiwan University} Taipei 10617, Taiwan, ROC Email: [email protected] Fax: 886-2-2363-3917.

al., 1998) and transversal configurations (Robertson and O’Connor, 1986; Herrmann and M¨uller, 1995) have been investigated in the literature, and rotation can be applied as an additional control means.

With an induction heated needle-eye configuration, some improvement of radial dopant distribution by rotation (up to 20 rpm) were reported under a weak magnetic field (0.1-0.14T) by Riemann et al. (1996). With a highly simplified model, Ma et al. (2000) further investigated the role of a weak transversal magnetic field and rotation on the flow structure. They claimed that the centrifugal pumping force overwhelms the inward convection and produces radial outward flow over the entire crystal-melt interface leading to an improvement of radial segrega-tion. However, for a small-scale FZ zone, the interac-tion of Marangoni flow and rotainterac-tion could be quite dif-ferent. Under magnetic fields the effect of rotation is not yet fully understood.

Self-consistent simulations must be regarded as a good means to shed some light on this problem and, in particu-lar, to provide a detailed picture of the transport behavior. In a recent study by Lan and Yeh (2004), the effects of axial and transversal magnetic fields were investi-gated numerically through a three-dimensional (3D) self-consistent model. They illustrated that the transver-sal field is more effective in suppressing the unsteady Marangoni flow. The required magnetic strength for get-ting a steady flow by means of the transversal field is lower than that required by the use of the axial one. However, the flow suppression/dumping by the transver-sal field tends to be effective only in the meridian plane parallel to the magnetic field. The flow in the plane per-pendicular to the direction of the magnetic field is not suppressed due to the induced electric potential, and this leads to a highly asymmetric molten zone and growth interface. The use of rotation, especially the counter-rotation, is believed to be an useful strategy for reduc-ing such an asymmetry. The growth experiment for the

needle-eye configuration (Riemann et al., 1996) is a pos-itive and relevant example.

To simulate the FZ growth problem in a self-consistent manner, because the zone shape is unknown and cou-pled with the heat flow, numerical methods able to pro-vide an ”ensemble” picture of the system are required. During the past 10 years, numerical methods based on a finite-volume method have been developed by our re-search group (Lan and Yeh, 2004; Lan, 1996; Lan and Liang, 1997). The global Newton’s method (Lan, 1996) has been found powerful for 2D problems. However, for 3D problems, it is less robust and requires tremen-dous computer memory and CPU time (Lan and Liang, 1997). Along these lines, point iterative methods seem to reduce the memory dramatically required by 3D mod-eling; however they suffer slow convergence; when the grid is refined the convergence rate degrades rapidly, be-cause long-wavelength errors cannot be removed effec-tively. The slow convergence associated with the solu-tion of the transport equasolu-tions can significantly affect the global convergence of the zone shape computation. How-ever a multigrid (MG) approach has been recently devel-oped (Lan and Yeh, 2004; Lan and Liang, 1999) that can amend this slow convergence problem.

Recently, Lappa (2003, 2004, 2005) and Gelfgat et al. (2005) also developed efficient numerical methods for solving 3D FZ problems. However, these studies did not take into account melt/solid front tracking, i.e. a fixed geometry was considered.

In this report, we investigate the effects of rotation for a small-scale FZ silicon growth under axial and transver-sal static magnetic fields through 3D self-consistent sim-ulation relying on the aforementioned efficient finite vol-ume/multigrid method (Lan and Liang, 1999). The flow structure, dopant distribution, and zone shape for various growth conditions are presented and discussed. The next section describes the mathematical model. Section 3 is devoted to the numerical method followed by the results and discussion. Conclusions are drawn in Section 5.

2 Mathematical model

A generic FZ crystal growth in a mono-ellipsoid mirror furnace is illustrated in Fig. 1; the total sample length is L. The mesh used for calculation (after convergence) is also shown. The long axis of the mirror (a) is 9 cm and the short axis (b) is 8 cm (Schweizer et al., 1999;

Figure 1 : Schematic of a floating-zone crystal growth

in a mono-ellipsoid mirror furnace; a sample mesh for simulation is also shown; a small portion of the melt is removed to illustrate the mesh inside the molten zone.

Lan and Yeh, 2004; Lan and Chian, 2001). Since ax-isymmetry is no longer assumed here, the system is de-scribed by a Cartesian coordinate (x,y,z) system. The point-source model (Lan, 1996) for the mirror furnace is adopted for the heat input condition. A more complete model for the lamp heating has been reported by Rivas and Haya (2002), but it requires tedious computation for an arbitrary zone shape.

The feed rod is steadily moved downward at speed Uf,

while the crystal pulling speed is Uc. For a pseudo-steady

state growth, for the same density of the feed and the crystal, AfUf = AcUc, where Af and Ac are the

cross-section area of the feeding rod and of the crystal, respec-tively. However, the shape of the grown crystal, which can be represented by a local radius Rc(θ), remains

un-known and needs to be calculated;θ is the azimuthal an-gle. When the molten zone is asymmetric, the cross sec-tion of the steadily-growing crystal is not exactly round.

Then, rotating the non-axisymmetric crystal can generate time-dependent oscillations. For the sake of the simplic-ity, we have ignored this effect; this is reasonable when the deviation from the round shape in the cross section is small.

The flow and temperature fields, as well as the melt/feed (hf(y,z)) and melt/crystal (hc(y,z)) interfaces, and the

meniscus (Rm(x,y,z)) are also represented in Cartesian

coordinates.

In this study, growth conditions have been chosen to have a steady-state growth, so that time-dependent calculation is not necessary.

The dimensionless variables are defined by scaling length with the feed rod diameter Df, velocity with αm/Df,

pressure with ρmα2m/D2f, temperature with the melting

point Tm, and dopant concentration by the feed rod

con-centration C0, where αm is the melt thermal

diffusiv-ity and ρmthe melt density. The steady-state governing

equations describing the convection, with the Boussinesq approximation, and heat and dopant transport in the melt read: ∇· vvv∗_{= 0, (1)} vvv∗· ∇vvv∗= −∇P∗+Pr∇2vvv∗+FFF∗, (2) vvv∗· ∇T∗= ∇2T∗, (3) vvv∗· ∇C∗=Pr Sc∇ 2_C∗_{, (4)}

where the body force term FFF∗in Eq. (2) (Baumgartl and M¨uller, 1992; Ben Hadid and Henry, 1996) can be writ-ten as follows:

FFF∗= PrRaT(T∗−1)eeex+PrHa2(−∇Φ∗+vvv∗×eeeB) ×eeeB

In the above equations, vvv∗, P∗, T∗, and C∗ are the di-mensionless velocity, pressure, temperature and dopant concentration, respectively; without the asterisk, the vari-ables are assumed to be dimensional. In addition, Pr= νm/αmis the Prandtl number and Sc= νm/D the Schmidt number, whereνmis the kinematic viscosity and D is the

dopant diffusivity. Also, eeex and eeeBare the unit vectors in

the x and magnetic directions, respectively, and they are considered in the x−y plane only:

eeeB= cos(β)eeex+sin(β)eeey,

where eeey is the unit vector in the y direction. In this

study,β = 0◦and 90◦are used for the axial and transver-sal fields, respectively.

The associated dimensionless thermal Rayleigh (RaT)

and Hartmann (Ha) numbers are defined as:

RaT =

βTg0TmD3f

νmαm

; Ha= |B|Df(σm/µm)1/2,

whereβTis the thermal expansion coefficient of the melt, g0 the gravitational acceleration, σm the melt electric

conductivity, µmthe melt viscosity, and| B | the applied

magnetic field strength. Because the solutal effect on the solidification temperature and the melt density has been neglected due to the light-doping approximation, the value of C0is arbitrary.

Furthermore, Φ is the electric potential, which requires the solution of the continuity equation for the electri-cal current density vector jjj∗, i.e., ∇ · jjj∗ = 0, where

jjj∗= σ∗_i(−∇Φ∗+ vvv∗× eeeB); σ∗_i = σi/σm, i= (m,c, f ), is the normalized electric conductivity of the melt (m), the crystal (c), and the feed rod ( f ). Here, the electrical cur-rent density vector jjj∗has been rescaled byσmαm|BBB|/Df

andΦ∗by αm|BBB| . In other words, we need to solve the following equation as well,

∇· σ∗

i∇Φ∗= ∇· σ∗i(vvv∗×eeeB). (5)

The magnetohydrodynamic (MHD) equations used above are the same as the MHD2 model in Baumgartl and M¨uller (1992) and the model used by Ben Hadid and Henry (1996). Furthermore, due to the light-doping ap-proximation, the interfaces are assumed isothermal and at the melting temperature corresponds to that of pure silicon. Accordingly, thermoelectromagnetic convection (TEMC) (Khine et al., 2000) can also be neglected. In the crystal (c) and the feed rod ( f ), only heat transfer needs to be considered:

∂T∗

∂τ +vvv∗i · ∇T∗= ∇· κi∇T∗,i = ( f ,c) (6)

where vvv∗_i is the dimensionless solid velocity including the rotational component. Also,κi= αi(T)/αmis the

dimen-sionless thermal diffusivity of feed and crystal;αi is the

thermal diffusivity of the feed rod(i = f ) or of the crystal (i= c) .

The no-slip condition is used for the velocity at solid boundaries. At the free surface, the shear stress balance is imposed:

where τ : nnnsss is the shear stress at the n − s plane of the free surface; n and s are the unit normal and tangential vectors at the free surface, respectively. Also, Ma is the Marangoni number, defined as:

Ma=|∂γ/∂T|TmDf

ρmvmαm

where ∂γ/∂T is the surface-tension-temperature coeffi-cient of the melt. Two tangential directions need to be considered for the stress balance. In addition, the kine-matic condition (nnn·vvv∗= 0) at the free surface and the

nor-mal stress balance (the Young-Laplace equation) must be also satisfied, i.e.,

τ∗_{: nnn nnn= (2H}∗_{)Bo +λ}∗

0, (8)

where 2H∗is the mean curvature scaled by 1/Df, Bo=

γ/(ρmg0D2f) is the static Bond number; γ is the surface

tension of the melt. The detailed procedure for calculat-ing the mean curvature can be found elsewhere (Lan and Liang, 1997). Also,λ∗₀ is a reference pressure head that needs to be determined to satisfy the growth angle con-straint for a steady growth, i.e., at the melt/gas/crystal tri-junction line,

nnnm· nnnc= cos(φ0), (9)

where nnnmand nnncare the unit vectors at the melt and

crys-tal surfaces, respectively, at the tri-junction line and φ0

is the growth angle for the growing crystal having a con-stant local radius. The growth angle can also depend on the crystallographic orientation. During calculation, the radius and the shape of the crystal are not known a

pri-ori. Therefore, the growth angle constraint is used for

the calculation of the local crystal radius Rc(θ).

How-ever, the static headλ0is used to satisfy the global mass

conservation, Ac= AfUf/Uc. With rotation, the interface

velocity vvv= (Uieeex+ rΩieee_θ)/γc, i= ( f ,c), where Uf and Uc are the feeding and pulling speeds, respectively; Ωf

and Ωc are rotations speeds of the feed rod and crystal

rod, respectively; γc is the density ratio of the crystal to

the melt.

The thermal boundary conditions at the growth and feed-ing fronts are set by the heat flux balances:

Q∗|m−Q∗|i+γcSt(vvv∗i · nnn) = 0,i = (c, f ) (10)

where nnn is the unit normal vector at the feeding or

grow-ing interface pointgrow-ing to the melt. Q∗|m, Q∗|c, and Q∗|f

are the dimensionless total heat fluxes at the melt, the crystal, and the feeding sides, respectively. The Stefan number St= ∆H/CpmTm scales the heat of fusion∆H released during solidification to the sensible heat in the melt; Cpmis the specific heat of the melt. The heat

trans-fer between the sample and the surrounding is similar to that in the 2D case described before (Lan, 1996). The boundary condition for the dopant transport at the growth interface is given by (where the solid-state diffusion is neglected):

n

nn· ∇C∗|m+ Sc

Prγc(1 −K)C

∗_(vvv∗_{· nnn) = 0, (11)}

where K is the segregation coefficient. For the feeding interface,(1 −K)C∗is replaced by(C∗−1).

The boundary condition for the electrical potentialΦ is set by the continuity of the electrical current density at the interfaces. In addition, the surface of the material and the end boundaries are assumed electrically insulated, i.e., jjj∗· nnn = 0, where nnn is the surface normal. Again, because jjj∗= σ∗_i(−∇Φ∗+vvv∗×eeeB), the insulation condi-tion becomes nnn· ∇Φ∗= nnn · (vvv∗×eeeB).

3 Numerical method

3.1 Co-ordinate transformation

Because the melt/solid interfaces and the meniscus are not know a priori, the co-ordinates (x,y,z) are trans-formed into a general curvilinear co-ordinate (ξ1,ξ2,ξ3) system, that fits all the interfaces and the free surface as shown in Fig. 2. Simple geometric transformations are used as follows: crystal: r= Rc(ξ2)ξ 3 (ξ3_{) (12)} θ = 2πξ2(ξ2_{) (13)} y= r cos(θ); z= r sin(θ) (14) x= hc(ξ2,ξ3)ξ 1 (ξ1_{) (15)} melt: r= Rm(ξ1,ξ2)ξ 3 (ξ3_{) (16)} x= hc(ξ2,ξ3) +[hf(ξ2,ξ3) −hc(ξ2,ξ3)]ξ 1 (ξ1_{) (17)}

Figure 2 : Schematic illustration of co-coordinate transformation: (a) the physical domain, and (b) the computational domain. feed rod: r= Rf(ξ2)ξ 3 (ξ3_{) (18)} x= hf(ξ2,ξ3) +[L−hf(ξ2,ξ3)]ξ 1 (ξ1_{) (19)}

whereξi(ξi) is a stretch function ranging from 0 to 1 in the computational domain. Simple analytical functions could be used to control grid distribution. For example, a hyperbolic tangent function with the following form for ξ1(ξ1_{) is used for the melt:}

ξ1(ξ1_{) = 0.5
1+tanh}B[(ξ 1_−ξ1 a)/(Nξ1−ξ 1 a) −0.5 0.5B  , (20) where B is a stretch constant and N_ξ1 the number of

con-trol volumes (CVs) in theξ1direction. During iterations, if the melt/ambient free surface shape Rm, the grown

crystal shape Rc, and the melt/solid interfaces (hf and hc) are estimated based on the proper boundary

condi-tions, the above coordinate transformations can be used to define the boundaries of CVs. For example, the coordi-nates of each CV in the physical space, as shown Fig. 3a,

are calculated according to the transformations. The cell faces are represented by surface vectors, i..e., A1, A2, and

A3, while the values of variables in each CV are defined at the geometric center.

3.2 Finite volume integration

After coordinate transformations, following Lan and Liang (1999) the governing equations can be rewritten in a general conservation-law form for each block in (ξ1, ξ2_,ξ3_{) :}

∂

∂ξi(Ciφ +Di) +JSφ= 0, (21)

where J= det|∂(x,y,z)/∂(ξ1,ξ2,ξ3)|. The generic vari-able φ stands for the Cartesian velocities (u∗, v∗, w∗ in the x-, y-, z-components, respectively), temperature and dopant concentration, respectively. For the equation of continuity, the variable φ can be set to unity. The de-tailed coefficients of Eq. (21) can be found in Lan and Liang (1999). Then, the finite volume method consists of simply integrating Eq. (21) over each CV in the compu-tational domain. After the Gauss theorem is applied, the

integration over each CV results in a flux balance equa-tion:

Ie−Iw+In−Is+It−Ib+ 

∆VSφdV = 0, (22)

where Ii, i= (e,w,n,s,t,b), represents the fluxes of φ

cross the faces of the CV. Each of the fluxes Iiis made of

two distinct parts, namely a convective contribution Ciφi

and a diffusive contribution Di. Both terms are

approxi-mated by the central difference scheme. In fact, accord-ing to the conditions considered here, upwind schemes are not necessary.

Owing to the fact that the pressure does not appear ex-plicitly in the equation of continuity, the use of linearly interpolated velocity from the face values for the col-located grids could lead to the pressure/velocity decou-pling, i.e., the so-called checkerboard pressure oscilla-tion. To amend this, we have adopted the Rhie-Chow mo-mentum interpolation scheme (Rhie and Chow, 1983) to interpolate the velocity values at cell faces directly from the momentum equations. Then, the SIMPLE method (Patankar, 1980) is used to correct the pressure for satis-fying the continuity equation.

3.3 Solution scheme

The finite volume approximation in Eq. (22) leads to a general nonlinear form at each CV center as

aφ_Pφ −

∑

aφ_nbφnb= Sφ∆V (23)

for variableφ, and φ = v∗,P∗,T∗ and C∗. The assembly of the above equation for all variables and all blocks leads to a sparse block matrix. The strong implicit procedure (SIP) based on the incomplete LU decomposition (Stone, 1968) is used to solve the block matrix. Then global iter-ation is performed variable wise and block wise. Never-theless, the performance is found satisfactory for coarse grids, but becomes poor as the grid is refined. In fact, for most of iterative methods, the convergence becomes slower as the number of equations increases. The con-vergence rate is usually fast at the initial stage, but it degrades rapidly as the iteration proceeds further. As mentioned before, this is because, generally, the short wavelength errors can be smoothed out efficiently by lo-cal iterations, but when fine grids are used, the reduc-tion/control of long wavelength errors becomes more dif-ficult.

An effective way to overcome the slow convergence is to use different meshes during iterations, i.e., the so-called multigrid (MG) methods. Because we are dealing with a free boundary problem, an accurately description of the boundary shapes is crucial. Therefore, the interface updates have to be considered for the finest grid. After setting up the finest grid, eight contiguous control umes (the kid cells) are coalesced into a larger finite vol-ume (the parent cell). A simple parent-kid link-list data structure is used for computation. With a given boundary shape, the full approximation scheme (FAS) proposed by Brandt (1977) is adopted here. After the MG iteration convergence process, the melt/solid interfaces are calcu-lated using the isotherm migration method (IMM). Then, the meniscus calculation is performed according to the normal stress balance, in which there is an outer iteration loop for moving the local radii and getting the reference pressure to satisfy the growth angle and the overall mass conservation. A conceptual sketch of this V-cycle itera-tion procedure is shown in Fig. 3b. Usually, the total in-ner iteration number for the fine grid, required to reduce the relative errors of about five orders, is of several hun-dred. This also includes about 30 to 120 iterations for the interfaces. One calculation takes about 2 hours of CPU time using a Pentium-4 (2GHz) PC. A sample converged mesh for calculation is shown in Fig. 1. As shown, finer grid spacing is placed near the interfaces to enhance the accuracy of calculation. Detailed benchmark comparison with the previous 2D and 3D calculations can be found elsewhere (Lan and Yeh, 2004; Lan and Liang, 1999).

4 Results and discussion

The growth of an 8-mm Si crystal with phosphorous dop-ing is considered; the heatdop-ing configuration is illustrated in Fig. 1 and the power is set at 700W (Lan and Yeh, 2004). Although the furnace model is highly simplified, the power consumption is still comparable to real ex-periments (500∼600W, Dold et al., 1998). The detailed physical properties and input conditions are the same as in Lan and Yeh (2004).

Fig. 4 shows the calculated results with 5 rpm counter rotation for an axial magnetic field at 0.5T. At this condi-tion, a steady-state result is obtained. The results on the

x−y and x−z planes show that the flow in the core region

is greatly suppressed, but the thermocapillary flow is still quite strong close to the melt/gas surface. This is consis-tent with previous 2D simulations (Lan, 1996). The flow

Figure 3 : (a) a typical control volume; (b) the schematic

of V-cycle multigrid iteration with the interface and free boundary updates at the top level.

and thermal fields on the x− z and x − y planes seem to be identical in terms of symmetry. However, if we exam-ine the result at x= 4.1cm on the y − z plane (about at the middle of the molten zone), the temperature exhibits a maximum at given locations; due to the related ther-mal gradients, there are eight flow vortices caused by the thermocapillary force. As a result, the induced electri-cal potential distribution on the same plane also exhibits an eight-cell structure. These results are very similar to those in the case of no rotation (Lan and Yeh, 2004), with the exception of the potential field.

When the magnetic field strength is increased to 1T, the result becomes perfectly axisymmetric. The electric po-tential contours on the plane at x= 4.1cm are perfect cir-cles having the maximum field at the center. Increasing the rotation speed up to 20 rpm does not change much the overall flow and thermal fields.

Figure 4 : Calculated flow and thermal fields at different

cross sections, as well as the electrical potential at the plane of z= 4.1cm (lower right figure), in the molten zone under an axial magnetic field of 0.5T with 5 rpm counter rotation.

With a transversal field, a steady-state result is obtained at a lower magnetic field strength (0.15 T). The calcu-lated result for 0.5T at a counter rotation of 5 rpm is illustrated in Fig. 5. As shown, the results are asym-metric, but have a two-fold symmetry viewed from the

y− z plane. The surface zone length on the plane (x − z

plane) perpendicular to the magnetic direction is longer than that on the parallel plane. As mentioned previously, the melt flow on the plane perpendicular to the magnetic field is not suppressed as that on the parallel plane (lead-ing to a longer zone length there). In the y− z plane (at

x= 4.1cm), the isotherms have an ellipsoid shape, while

the maximum temperature appears at the surface of the

Figure 5 : Calculated flow and thermal fields at different

cross sections, as well as the electrical potential at the plane of z = 4.1cm (lower right figure), in the molten zone under an transversal magnetic field of 5T with 5 rpm counter rotation.

of Lan and Yeh (2004), it is clear that 5 rpm counter ro-tation has little effect on the overall flow and the thermal fields. The calculated potential field is also similar to that obtained in the case of no rotation.

Interestingly, although the effect of rotation on the ap-parent heat flow for both axial and transversal fields is weak, its effect on the zone shape and dopant segrega-tion is not trivial. Figure 6 shows the effect of rotasegrega-tion on the asymmetry of the zone length at the surface; the asymmetry degree is herein defined in terms of the max-imum difference of the zone length at the melt surface. Under an axial field, the heat flow is axisymmetric (or nearly), and the zone length is uniform. However, under a transversal field, the zone length is not uniform

with-Figure 6 : Effect of rotation on the asymmetry of the

zone length for axial and transversal fields.

out rotation, and the asymmetry degree increases with the transversal field magnitude. However, the asymme-try reduces rapidly with rotation. As the rotation speed is greater than 5 RPM, the zone length becomes uniform. The reason for this can be found in the Stefan boundary condition in Eq. (10). With rotation, the Stefan term in the angular direction, i.e., St(r Ωf,ceeeθ· nnn), plays an

im-portant role in reducing the asymmetry of the interface shapes. Similar situations were also discussed in hori-zontal zone melting (Lan et al., 2000a) and Bridgman crystal growth (Lan et al., 2000b).

The effect of rotation on the radial dopant segregation in terms of the normalized radial concentration difference is summarized in Fig. 7a. Without rotation, the case with 0.5T transversal field has the lowest radial segrega-tion. As also discussed in earlier analyses (Lan and Yeh, 2004), this is due to better dopant mixing caused by the stronger Marangoni convection in the plane perpendicu-lar to the magnetic field (that is not suppressed). Interest-ingly, the radial segregation increases with rotation when the speed is less than 5 rpm. As the rotation speed is further increased, the radial segregation decreases again. Surprisingly, the case under 1T transversal field shows a different trend. However, for the axial field, the effect of rotation is trivial; the flow at 1T is greatly dumped, and the flow due to rotation tends to be also suppressed. Clearly, the radial segregation is affected by the local mixing in the melt core near the growth interface, but is

Figure 7 : Effect of rotation on (a) radial dopant segregation and (b) effective segregation coefficient (Keff) for axial

and transversal fields.

difficult to discern these effects on the basis of the global flow structure shown in Figs. 4 and 5. Nevertheless, it should be pointed out that the poorer non-uniform dopant mixing is, the larger radial segregation becomes. A uni-form flow gives no mixing and thus the least segregation. Similar results can be illustrated and discussed on the ba-sis of the effective segregation coefficient Ke f f, which is

defined by C0/C (Lan et al., 2000b); C is the average

dopant concentration in the molten zone. A better mixed zone leads to a smaller Ke f f value. As shown in Fig. 7b,

for the rotation up to 20 RPM, in general, rotation has little effect on the dopant mixing. This is consistent with our observation. However, for the transversal field, the effect of rotation on the global mixing is not trivial at low rotation speeds. This also explains the subtle change of the radial segregation caused by slow rotation shown in Fig. 7a.

Finally, the 3D view of the dopant contours shown in Fig. 8 allows to easily discern the related effect of rotation. For the transversal fields, the rotation brings the dopant toward the center of the interface bridging the two dopant peaks. On the other hand, similar to the Stefan effect at the interface shape, rotation also generates angular segre-gation due to the asymmetric interface shape, particularly near the crystal surface, where the dopant concentration has the lowest value. As a result, the rotation tends to

increase the radial segregation for the transversal field. Similar simulation results for Bridgman crystal growth could be found elsewhere. For the axial field at 1T, the effect of rotation is not significant.

For silicon, the electric conductivity in the crystal is about one order smaller than that in the melt. Therefore, the centrifugal melt pumping (Hjellming and Walker, 1986; Lan and Yeh, 2004b) due to rotation in the transversal field is not significant in the present study. However, when the dopant level is high, the solid electri-cal conductivity can be increased significantly. There-fore, we also performed calculations by letting σc =

σm; accordingly the axial electric current across the

melt/crystal interface is increased significantly and en-hances the centrifugal melt pumping. As shown in Fig. 8, such an enhanced centrifugal pumping induces an over-shoot on the dopant field at the center of the interface. As a result, the radial segregation increases.

5 Conclusions

In this study, we have performed 3D calculations to in-vestigate the effect of rotation on the flow and dopant fields, as well as the zone shape, for a small-scale floating-zone silicon growth in a mono-ellipsoid mirror furnace under axial and transversal magnetic fields. In general, when the counter rotation speed is less than 20

Figure 8 : Effect of rotation on dopant distribution at the growth interface for transversal and axial magnetic fields;

the height and the color of each plot indicate the dopant concentration. The upper two rows are for the transversal field, while the lowest row is for the axial field. The cases havingσf = σc= σmare also included for comparison.

rpm, the rotation effect on the overall flow and thermal fields is small. However, the asymmetry of the zone length and growth interface due to the transversal field can be significantly reduced by rotation. The effect of rotation on the radial segregation is also significant for the transversal field, while its effect is smaller for the ax-ial field case. Increasing the electrical conductivity of the crystal enhances centrifugal melt pumping. For the transversal field at a slow rotation rate, the centrifugal pumping significantly increases radial segregation. Nev-ertheless, rotation has little effect on the global dopant mixing.

References

Amberg G. and Shiomi J. (2005): Thermocapillary flow

and phase change in some widespread materials pro-cesses, FDMP, vol. 1, no. 1.

Baumgartl, J.; M ¨uller, G. (1992): Calculation of the

effect of magnetic field damping on fluid flow, - compar-ison of magnetohydrodynamic models of different com-plexity, In Proc. 8th Symposium on Materials and Fluid Science in Microgravity (No. SP 333.), ESA Publication,

Brussels, pp. 161-165.

Ben Hadid, Henry, D. (1996): Numerical simulation of

convective three-dimensional flows in a horizontal cylin-der uncylin-der the action of a constant magnetic field, J.

Crys-tal Growth, vol. 166, pp. 436-445.

Brandt, A. (1977): Multi-level adaptive solutions to

boundary value problems, Math. Comp., vol. 31, pp. 333-390.

Dold, P.; Cr¨oll, A.; Benz, K. W. (1998): Floating-zone

growth of silicon in magnetic fields. I. Weak static axial fields, J. Crystal Growth, vol. 183, no. 4, pp. 545-553.

Gelfgat A. Yu., Rubinov A., Bar-Yoseph P. Z. and Solan A., (2005): On the three-dimensional instability of

thermocapillary convection in arbitrarily heated floating zones in microgravity environment, FDMP, vol.1, no. 1.

Herrmann, F. M.; M ¨uller, G. (1995): Growth of 20 mm

diameter GaAs crystals by the floating-zone technique with controlled As-vapour pressure under microgravity,

J. Crystal Growth, vol. 156, pp. 350-360.

Hjellming, L. N.; Walker, J. S. (1986): Melt motion

in a Czochralski crystal puller with an axial magnetic field, isothermal motion, J. Fluid Mechanics, vol. 164, pp. 237-273.

Khine, Y. Y.; Walker, J. S.; Szofran, R. F. (2000):

Thermoelectric magnetohydrodynamic flow during crys-tal growth with a moderate or weak magnetic field, J.

Crystal Growth, vol. 212, pp. 584-596.

Lan, C. W. (1996): Effect of axisymmetric magnetic

fields on radial dopant segregation of floating-zone sil-icon growth in a minor furnace, J. Crystal Growth, vol. 169, no. 2, pp. 269-278.

Lan, C. W.; Liang, M. C. (1997): A three-dimensional finite-volume/Newton method for thermal-capillary problems, Int. J. Numerical Methods in

Engi-neering, vol. 40, pp. 621-636.

Lan, C. W.; Liang, M. C. (1999): Multigrid methods

for incompressible heat flow problems with an unknown interface, J. Comput. Phys., vol. 152, pp. 55-77.

Lan, C. W.; Chian, J. H. (2001): Three-dimensional

simulation of Marangoni flow and interfaces in floating-zone silicon crystal growth, J. Crystal Growth, vol. 230, pp. 172-180.

Lan C. W.; Chian, J. H.; Wang, T. Y. (2000a):

Inter-face control mechanisms in horizontal zone-melting with slow rotation. J. Crystal Growth, vol. 218, pp. 115-124.

Lan, C. W.; Liang, M. C.; Chian, J. H. (2000b):

In-fluence of ampoule rotation on three-dimensional con-vection and segregation in Bridgman crystal growth un-der imperfect growth conditions, J. Crystal Growth, vol. 212, pp. 340-351.

Lan, C. W.; Yeh, B. C. (2004): Three-dimensional

simulation of heat flow, segregation, and zone shape in floating-zone silicon growth under axial and transversal magnetic fields, J. Crystal Growth, vol. 262, pp. 59-71.

Lan, C. W.; Yeh, B. C. (2004b): Three-dimensional

analysis of flow and segregation in vertical bridgman crystal growth under a transversal magnetic field with ampoule rotation, J. Crystal Growth, vol. 266, pp. 200-206.

Lappa, M. (2003): Three-dimensional numerical

simu-lation of Marangoni flow instabilities in floating zones laterally heated by an equatorial ring, Physics of Fluids, vol. 15, no. 3, pp. 776-789.

Lappa, M. (2004): Combined effect of volume and

gravity on the three-dimensional flow instability in non-cylindrical floating zones heated by an equatorial ring,

Physics of Fluids, vol. 16, no. 2, pp. 331-343.

Lappa, M. (2005): Analysis of flow instabilities in

con-vex and concave floating zones heated by an equatorial ring under microgravity conditions”, Computers &

Flu-ids, vol. 34, no. 6, pp. 743-770

Ma, N.; Walker, J. S.; Ludge, A.; Riemann, H. (2000):

Silicon float zone process with a weak transverse mag-netic field and crystal rotation, J. Electrochemical

Soci-ety, vol. 147, pp. 3529-3534.

Nakamura, S.; Hibiya, T.; Kakimoto, K.; Imaishi, N.; Nishizawa, S.; Hirata, A.; Mukai, K.; Yoda, S.; Morita, T. S. (1998): Temperature fluctuations of the

Marangoni flow in a liquid bridge of molten silicon un-der microgravity on board the TR-IA-4 rocket, J. Crystal

Growth, vol. 186, no. 1-2, pp. 85-94.

Rhie, C. M.; Chow, W. L. (1983): Numerical study of

the turbulent flow past an airfoil with trailing edge sepa-ration, AIAA J., vol. 21, no.11, pp. 1525-1532.

S.V. Patankar (1980): Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington DC.

Riemann H., Ludge A., Hallmann B., Turschner T.,

High Purity Silicon IV, C.L. Claeys, P. Rai-Choudhury,

P. Stallhofer, J. E. Maurits, Editors, PV 86-13, The Electrochemical Society Proceedings Series, Pennington,

NJ, pp. 49 (1996).

Rivas, D.; Haya, R. (2002): Analysis of secondary

ra-diation (multiple reflections) in monoellipsoidal mirror furnaces, J. Crystal Growth, vol. 241, pp. 249-260.

Robertson, G. D.; O’Connor, D. J. (1986): Magnetic

field effects on floating-zone Si crystal growth. II. strong transversal fields, J. Crystal Growth, vol. 76, pp. 100-110.

Schweizer, M.; Cr¨oll, A.; Dold, P.; Kaiser, T. H.; Lichtensteiger, M.; Benz, K. W. (1999): Measurement

of temperature fluctuations and microscopic growth rates in a silicon floating zone under microgravity, J. Crystal

Growth, vol. 203, no. 4, pp. 500-510.

Stone, H. L. (1968): Iterative solution of implicit approximation of multidimensional partial differential equations, SIAM J. Numer. Anal., vol. 5, pp. 530-558.