Amplitudes of doping striations: comparison of numerical calculations and analytical approaches

j. . . . . . . .

ELSEVIER

ClRYSTAL OIROWTH

Journal of Crystal Growth 171 (1997) 373-379

Amplitudes of doping striations: comparison of numerical calculations and analytical approaches T. Jung, G. Miiller * lnstitut ffir Werkstoffwissenschaften LS VI, Uniuersit~t Erlangen-Nfirnberg, Martensstrasse 7, D-91058 Erlangen, Germany Received 23 May 1996; accepted 26 July 1996

Abstract

Transient, axisymmetric numerical calculations of the heat and species transport including convection were performed for a simplified vertical gradient freeze (Bridgman) process with bottom seeding for GaAs. Periodical oscillations were superimposed onto the transient heater temperature profile. The amplitudes of the resulting oscillations of the growth rate and the dopant concentration (striations) in the growing crystals are compared with the predictions of analytical models.

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

An important advantage of the dynamical vertical gradient freeze (DVGF) for the growth of semiconductors is the stability of the temperature profile. However, if temperature oscillations in the heating zones would occur, they induce oscillations of the growth rate and of the resulting dopant distribution in the grown crystal. Predictions of the amplitude of these dopant striations as a function of amplitude and frequency of disturbations in the heating profile, or the calculations of tolerable values, have to take into account the heat transfer between heater and crucible, possible effects of convective heat transfer, the release of latent heat at the growth interface and the depth to which the disturbance penetrates into the melt, which is varying with the frequency. Up to now, to our knowledge there is no work in the

* Corresponding author. Fax: +49 9131 858495; [email protected].

E-mail:

literature treating this problem quantitatively. Wilson [1,2] published some numerical results for the Czochralski growth configuration, but the variation of the growth rate was prescribed. In the works of Van Run [3], Favier and Wilson [4] and Hurle et al. [5] a temperature oscillation at the edge of the thermal boundary layer is prescribed, which is difficult to relate to an applied outer temperature profile. Thevenard et al. [6] published an analytical model, where the oscillation of the growth rate is related to the oscillation of the temperature gradient at the growth interface, which in practical cases is also difficult to estimate. Concerning the effect of a given oscillation of the growth rate on dopant incorporation into the growing crystal, mainly the analytical models of Thevenard et al. [6] and Garandet [7] have to be mentioned. In this paper, we want to show some results from transient numerical calculations of the DVGF growth of GaAs, where temperature oscillations of varying frequencies have been superimposed onto the heating temperature profile. The re-

0022-0248/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S 0 0 2 2 - 0 2 4 8 ( 9 6 ) 0 0 7 0 4 - X

T. Jung, G. Miiller / Journal of C~stal Growth 171 (1997)373-379

374

suiting oscillations of dopant incorporation into the crystal are compared to the predictions of the analytical models [6,7].

2.

N u m e r i c a l

model

The numerical model used for the calculations presented here has been described already [ 8 - 1 1 ] . The model is transient and axisymmetric, the Boussinesq approximation is applied, the growth interface is determined by a front-tracking method neglecting capillary effects and the change of the melting point with composition. Calculations have been performed using the following conditions: height of crystal and melt together = 10 cm, radius of melt/crystal = 2.5 cm; pBN crucible, wall thickness = 1 ram; anisotropic heat conductivity A: in axial direction 62.7 W / ( K - m), in radial direction 3.37 W / ( K . m); temperature gradient, 20 K / c m along the crystal, 0 K / c m along the melt; this profile is coupled with the outer crucible wall by a heat transfer coefficient of 300 W / K . m 2 (see Fig. 1). growth starts at a crystal length of 3 cm, the temperature profile is translated with a constant velocity of 1 / z m / s for a time of 20000 s; after 20 000 s (2 cm of growth), a temperature fluctuation of the form AT = 0.5 sin[(27r/~')t] is

Table 1 List of the material parameters used for GaAs Material parameter

Value

Thermal expansion coefficient, ot Thermal conductivity (liquid), A: Thermal conductivity (solid), A~ Heat capacity (liquid and solid), Cp Kinematic viscosity, v Density, p Latent heat Distribution coefficient, k Diffusion coefficient dopant, D Prandtl number, Pr

1 × l0 -4 1,/K

17.8 W / ( K . m) 7.12 W / ( K . m) 434 J / ( K . kg) 4.88× 10 -7 m 2 / s 5.72× 103 k g / m 3 97 k J / m o l 0.1 I × 10 - s m 2 / s 0.068

added to the temperature profile with ~- varying from 1 to 5000 s; additionally, calculations are performed taking

t i• [ ......

temperalure[KI

:i. " . . . . . . B

1 t ~ ~"• ',',', ~,~ 1/', IX t1 l~ • ~ |/~11 ! ~ /~ ~ ' ' i i ~

-~ /~ i;~

l~X,

s s ,' : i j

,

.. ,

.• /

:

~ li~/, , f . / . : I r

i

1430

E 146(I Ct F1144014911 1470 G 1480

i ~s<,ol K 15101

L 1521)I M I530 I

cl~llCeulfrafitnlc/cr~

1540 1520

?

healer temperalure _ g_nitmle_ /

p_oi_nL--151_1K ~ ~-///__Z~/_ _

.............

~'1500 1480

~ 1460 1440 1420 1400 0.0

0.02

0.04 0.06 axial pos. [m]

0.08

0.1

Fig. 1. Temperature profile at the heater and the resulting temperature profile at the inner crucible wall.

Fig. 2. Temperature field, flow velocities and isoconcentration lines in the melt after 20000 s of growth, before applying oscillatory boundary conditions• Only one of four calculated velocity arrows has been drawn in axial resp. radial direction.

T. Jung, G. Miiller / Journal of Cr3,stal Growth 171 (1997) 373-379 O.3 [

/

A t I I I ,," : ~ •

at axis o f synuneuy . . . . ._ ._. . . . . . . . . . . ......

0.25 /., (I.15

<•1o -2

"~

:: ' i ,, q

vVVVl

Ix: .:." " (I

5(XI0

I(I(ICX)

15(X)0

20(X)0

[

O I

melt grid: 20x40

-

_~

t"

2

~

0.(15 (I.(I

375

l(i-~

2500{)

{~

5

time [sl

I0

15

2(J

25

radial position [mini

F i g . 3. Growth rate near the crucible wall and at the axis of

symmetry for a disturbance period of r = 1 0 0 0 s.

into account only conduction of heat to get an idea of the influence of convection. For material parameters, see Table 1.

5. Radial dependance of the amplitude of striations for numerical grids with 2 0 × 4 0 and 4 0 × 8 0 control volumes in the melt, oscillation period r = 1 0 0 s.

Fig.

the other disturbance periods, only the amplitudes of growth rates and dopant striations will be shown as functions of the radial position.

3.1. Influence of grid spacing on the numerical solution

3. Results of numerical calculations Fig. 2 shows temperatures, flow velocities and the dopant distribution after 20 000 s of growth. At this time the temperature disturbances are switched on. The Grashof number, calculated from the radial temperature difference in the melt at the height where the melt-solid interface is in contact with the crucible is 6.7 × 10 4. The Reynolds number, based on the maximum fluid velocity and the crucible radius is 40. Figs. 3 and 4 show, as an example, the growth rate and axial dopant distribution in the crystal for a period of temperature oscillations of r = 1000 s. For

The results presented here were obtained with a numerical grid of 40 × 80, 40 X 40 and 10 X 120 control volumes for the melt, crystal and crucible, respectively. For a strict check for the dependence of the solution on the numerical grid it would be necessary to repeat the calculations on a finer grid. Unfortunately, this is not possible because calculations take too much time then. Nevertheless, Fig. 5 shows that results for the amplitude of the striations in a selected example case are at least similar on the coarser grid with 20 X 40, 20 X 20 and 5 × 60 control volumes. Calculations with an even coarser grid

I).26 .... =.....'...."

"B O.24

@p4

al axis o f symmetry

=© 0.22 O.2

.........................

.~

,=- 0.18 '~

/,.,'

0.16 0.]4

o IOx20control volumes ill the melt 20x40

;> <:l i0 ~ 5 2 ill -'~

eaa 2

-

a~ crucible wall

ca "~ Ill 3 .-=_ s

-

Ikl2 0.1 5000

I(XI(IO

15000

2000~1

25(1(R1

time [s]

10-4 0

5

Ill

15

20

25

radial position [mm]

F i g . 4. Dopant distribution in the grown crystal near the crucible

F i g . 6. Radial dependence of the amplitude of growth rate for

w a l l and at the axis of symmetry for a disturbance period of

numerical grids with lO × 2 0 , 2 0 × 4 0 a n d 4 0 X 8 0 control volumes in the melt, 7 = lO s.

r=

1 0 0 0 s.

T. Jung, G. Miiller/Journal of Crystal Growth 171 (1997) 373-379

376

and taking into account convection and segregation are not possible. The convective roll at the phase boundary is then no longer resolved. However, they can be done for the temperature field alone, which could be of interest for a fast analysis of the response of the growth rate to varying disturbances of the heating profile. As an example, Fig. 6 shows the amplitude of the growth rate for a period of ~-= 10 s calculated with 10 X 20, 20 X 40 and 40 X 80 con-

T:lOs

trol volumes in the melt, without convection. The coincidence on all grids is excellent.

3.2. Response of the conuectiue flow to the temperature distribution To visualize the changes in the flow velocities induced by the oscillating boundary conditions, we have plotted in Fig. 7 streamlines calculated from the

T= 100s

T= I O00s

s t r e a m lines: a b c d e f g h i k I m

-10" -10" -10' -10" -10 ~ - 1 0 .... +10 '° +10" +10 ~ +10" +10' +10 +

o_. 3 3

~

~

J

Y

Fig. 7. Streamlines calculated from differential velocity changes with time for oscillation periods of 10, 100 and 1000 s. The dot on the sine curve on the left represents the phase of the applied temperature oscillation.

7". ,lung, G. Miiller / Journal of CO,stal Growth 171 11997) 373-379

differential changes of flow velocity with time. The arrows at the streamlines show the direction of the flow changes. Fig. 7 shows the appearance of inertial effects: For r = 1000 s, the changes in the flow velocity are in phase with the outer temperature oscillation (increasing in a clockwise direction while the outer temperature increases and vice versa). For r = 10 s, there is a phase shift of about - 7r/2 to - ~-, while the streamlines for r = 100 s show a more complicated behavior, indicating some sort of transition. As the sign of the radial temperature gradient does not change during an oscillation period, this must be an effect of inertia. Of course, one cannot expect a real melt to behave exactly the same way, there will certainly be three-dimensional effects. A more detailed discussion of this phenomenon would be beyond the scope of this paper, but this point had to be mentioned at least because the transition will be reflected in the results for the amplitudes of growth rate and striations presented in the following.

3.3. Amplitude of growth rate." numerical results To check for the influence of convective heat transfer, calculations with disturbance periods of 1, 10, 100 and 1000 s have been performed with and without convection (see Fig. 8). There is a remarkable difference in the oscillation of the growth rate with and without convection at a period of 10 s, which does not appear for the other periods calculated here, and is probably related to the resonance effect described in Section 3.2.

,z

5

0.14

': 0.12 > ,j 0.1

a[ i-=1

.,

-E 0.0~

0.04

"-- I).1)2

,.-

,IL'

m""

5

HIe

..4 .....................

,o

,."

.,'"

•J"

~, 006 -

"'"'"'

.~7cm

.-"

o...'" at axis of symmetry

..

.:~.$m....:t" ....

(111 I111

2

2

s

l0 ~

2

disturbation period m Is] Fig. 9. Amplitude of growth rate as a function of the period of the temperature disturbance (symbols). The lines represent polynomial fits which in Section 3.5 are input into the analytical models for the dopant incorporation.

Except the region at r ~ 10 s, the differences between the amplitudes calculated with and without convective flow are not significant. Therefore, for a future analysis of the growth rate amplitude with varying boundary conditions and different forms of temperature disturbances, it should be possible to get an overview with much faster calculations neglecting convective flow. Such calculations can give reasonable results even on quite coarse grids, decreasing computation times again significantly, as shown in Fig. 6, for the case of r = 10 s with 10 × 20, 20 × 40 and 40 × 80 control volumes in the melt region. Fig. 9 shows the amplitudes of the growth rate, calculated with convection, at the axis of symmetry and at the crucible wall. These results will be used as input for the analytical models for the dopant amplitudes in Section 3.5.

"F=l(10/ls

: ...................................

"-<>'/.I-

"-;;= -'-": ;.-:.:;-k;;~;;: - --..-7-.]%-.- .-~,.ZT~7.T............ T = l O l l s

.=

377

e ................................

j

f,.i'I

/

z

{,°;i N

1,, 2 ~_. 1040

.......... ........................

5

2

without COllVeCtio[i ............ '-...... __ l/ - - - with conveclion ,

,

10

15

, a/

2(1

25

radial posilion [ram] Fig. 8. Radial dependencies of amplitudes of the growth rate for various periods r of oscillations, with and without convective heat transfer.

1(i 5 5

10 15 radial IX~Sin mm

20

2

25

Fig. 10. Amplitudes of dopant incorporation into the growing crystal as a function of the radial position in the crystal, for different periods r of the temperature disturbance.

378

T. Jung, G. MMler / Journal of Crystal Growth 171 (1997) 373-379 8c = 10nmn

3.4. Fluctuation of dopant concentration." numerical results

1 05- ~ at r=1.67c m ~ = l m m

Fig. 10 shows the calculated amplitudes of the fluctuations of the dopant concentration in the grown crystal as a function of the radial position. The minimum at r = 17 mm and a period r = 100 s is related to the transition of the response of the convective flow to the temperature disturbance described in Section 3.2.

~" .~

C~

VIl

7,0 1-(iA°/2Fq)[1-(1

+4iFq) 1/2]

)< ½[I + ( 1 +4iFq) 1/2] - ( l - k )

A° ' (1)

with A° = 6~VG/D and Fq = o 9 6 2 / D . The result of Thevenard et al. is [6] v Cd = ~ Vl I ,

(2)

with (1 - k ) mC I

1 + (V1°)2) sinh(A6c) 2i ~oD

D

5

I---a~evenurd I ~ numerical results I .

.

.

.

2

.

.

• .....

8,.=10mm

/ ~ .....

10 4

~

~ l m m

1o-~_

3.5. Comparison of numerical results with the predictions of the analytical models The final result of the analytical model of Garandet presented in Ref. [7] is the following equation for the relative amplitude of the fluctuation of the doping concentration:

; /7 1()--s //~

/'

~

-

,o- t ¢ / 5 1()t 2

5 10 2 _'

5 10 3 2

5 1(14

oscillation peri{~l "r Is] Fig. 11. Comparison of striation amplitudes from numerical calculations with the predictions of the analytical models for the center (below) and a radial position of r = 1.67 cm.

extent of the solute boundary layer, 6c. 6c can be estimated from the model presented by Camel and Favier (e.g. Ref. [12]), which yields values for 6c of some millimeters for a Grashof number of 105 (see Ref. [11]). We have used in this paper the amplitudes of the growth rate shown in Fig. 9). Fig. 11 presents the comparison of the numerically obtained striation amplitudes with the results of Eqs. (1) and (2). Curves are shown for the radial positions r = 0 and r = 1.67 cm and assuming boundary layer extents of 0.1, 1 and 10 mm. Both analytical models give practically identical results. At the axis of symmetry, 1.9 c~

"~-gl° [ c o s h ( / ~ c )

-

~1.8 "O

exp(Vi°rJ2D)]

41.7

lOJ

v' 0

]-~

.3

mmetry

1.6 1.5

Acosh(ASo) - -f-~ (1 - 2k)sinh(Ar~)

°l4

A=

¢i10) 12D

+-D-'

(3)

~'1.3 1.2

(4)

Unknown in Eqs. (1) and (2) are usually the amplitude of the oscillation of the growth rate, V1, and the

0

10 20 30 40 cfistance from growth interface [mm]

50

Fig. 12. Axial solute distribution in the melt (without imposed temperature oscillations) at the axis of symmetry and at r = 1.67 cm.

T. Jung, G. Miiller/Journal of Crystal Growth 171 (1997) 373-379

the numerical results are fit better assuming a thickness of the boundary layer of 10 mm, at r = 1.67 cm with a value of 1 mm. This agrees well with the calculated axial solute distribution in the melt shown in Fig. 12. 3.6. Conclusions

Generally, the coincidence between the results from the analytical models [7,6] and from the numerical simulations is very good, which justifies the simplifying assumptions necessary to derive them. A somehow more complicated behavior of the striation amplitudes has to be expected in the frequency range between 0.01 and 0.1 Hz, where inertial (or resonance) effects in the response to the convective flow to the outer temperature disturbances can lead to effects not predicted by the analytical models. For a practical use of the analytical models, it is necessary to know the amplitude of a growth rate oscillation as a function of varying boundary conditions and types of disturbances, for which no theoretical models exist. Fortunately, it appears to be possible to obtain these data from numerical calculations on coarse grids and neglecting convection, except at disturbance periods around 0.I Hz.

C~~ dopant concentration in the crystal without oscillation V~° undisturbed growth rate V~~ amplitude of growth rate k equilibrium distribution coefficient m slope of liquidus curve D diffusion coefficient of dopant in the melt r oscillation period

References [I] [2] [3] [4] [5] [6] [7] [8] [9]

4. Names and symbols 6c C~

solute boundary layer thickness amplitude of the fluctuation of the dopant concentration in the crystal

379

[10] [1 I] [12]

L.O. Wilson, J. Crystal Growth 48 (1980) 435. L.O. Wilson, J. Crystal Growth 48 (1980) 451. A.M.J.G. van Run, J. Crystal Growth 47 (1979) 680. J.J. Favier and L.O. Wilson, J. Crystal Growth 58 (1982) 103. D.T.J. Hurle, E. Jakeman and E.R. Pike, J. Crystal Growth 3 / 4 (1968) 633. D. Thevenard, A. Rouzaud, J. Comera and J.J. Favier, J. Crystal Growth 108 (1991) 572. J.P. Garandet, J. Crystal Growth 131 (1993) 431. D. Hofman, T. Jung and G. Miiller, J. Crystal Growth 128 (1993) 213. D. Hofmann, T. Jung, K. Sch~ifer, D. Zemke and C. Miiller, Mater. Sci. Eng. A 173 (1993) 67. T. Jung, PhD Thesis, University Erlangen-Niirnberg, 1993. T. Jung and G. Miiller, J. Crystal Growth 165 (1996) 463. D. Camel and J.J. Favier, J. Physique 47 (1986) 1001.