- Email: [email protected]
Thermal stress analysis of silicon bulk single crystal during Czochralski growth
Journal of Crystal Growth 125 (1992) 102—111 North-Holland
j
~
CRYSTAL
GROW T H
Thermal stress analysis of silicon bulk single crystal during Czochralski growth N. Miyazaki, H. Uchida, T. Munakata Department of Chemical Engineering, Faculty of Engineering, Kyushu UniLersity. 6-10-1 Hakozak,, Higashi-ku, Fukuoka-shi, Fukuoka-ken 812, Japan
K. Fujioka Mechanical Engineering Research Laboratory, Hitachi, Ltd. 502 Kandatsu-machi, Tsuchiura-shi, Ibaraki-ken 300, Japan
and Y. Sugino Kofu Works, Hitachi, Ltd., Nishiyahata, Ryuoh-cho, Nakakoma-gun, Yamanashi-ken 400-01, Japan Received 15 April 1992; manuscript received in final form 21 June 1992
The thermal stress analysis of a silicon bulk single crystal with a diameter of 6 or 8 inches is performed in the cases of the [001] and [1111pulling directions by using a three-dimensional finite element program developed for calculating thermal stress in a bulk single crystal during the Czochralski growth. Elastic anisotropy and temperature dependence of material properties are taken into account in this program. The temperature distribution and shape of a silicon bulk single crystal which are required for the thermal stress analysis are obtained from a computer program for a transient heat conduction analysis which is specialized for the Czochralski growth. The stress components obtained from the thermal stress analysis are converted into the parameters related with dislocation density. The time variations of these parameters are shown in this paper. The relation between these parameters and the shape of the crystal—melt interface is discussed.
1. Introduction The dislocations formed in bulk single crystals during the Czochralski (CZ) growth have adverse effects on the performance of electronic devices, It is known that the thermal stress during the CZ growth plays an important role in the formation of dislocations, Although a lot of studies [1—8] have been performed concerning the thermal stress analysis of bulk single crystals during the CZ growth, there have been few studies which deal with the thermal stress analysis of such large-sized bulk single crystals as industrial use during a long period of the CZ growth. 0022-0248/92/$05.00 © 1992
—
Miyazaki et a!. [9,10] have developed a threedimensional finite element program for the thermal stress analysis of single crystals, which considers the elastic anisotropy and the temperature dependence of the material properties. On the other hand, Nakayama and Masaki [11] have developed a computer program for a transient heat conduction analysis which can simulate the CZ growth process of a bulk single crystal by producing finite elements automatically and by considering heat radiation between the bulk single crystal and its neighboring environment, i.e., melt, heater and radiation shield. In the present study, the thermal stress analysis of a silicon bulk single crystal with a diameter
Elsevier Science Publishers By. All rights reserved
N. Miyazaki et al.
/ Thermal stress analysis
of 6 or 8 inches during a long period of the CZ growth is performed by using both the abovementioned computer programs for a thermal stress analysis and a transient heat conduction analysis. The stress components obtained from the thermal stress analysis are converted into the parameters related with dislocation density, i.e., dislocation density parameters. The time variations of these parameters are shown and the relation between these parameters and the shape of the crystal—melt interface is discussed.
of Si bulk single crystal during CZ growth
The resolved shear stress can be calculated by the transformation of stresses from the Cartesian coordinates used in a thermal stress analysis to a slip system of a single crystal. When 0~ex is equal to zero, no dislocation is expected to be generated. Speaking exactly, 0~tot is not a dislocation density parameter but a measure for a thermal stress effective to glide strain. In addition to o~ and ~ we also used another dislocation density parameter °‘mex [13], which is defined as (3) When °~mexis less than unity, no dislocation is expected to be generated. The dislocation density parameter 0~ex indicates whether or not the thermal stress level in a bulk single crystal exceeds 0~cRs’but does not indicate how many times the thermal stress level is as large as ~CRS’ Therefore, we used another dislocation density parameter 0’mex as a quantitative measure indicating the difference between the thermal stress level and ~
°mex
2. Method of analysis 2.1. Elastic anisotropy In the present analysis, the elastic anisotropy is taken into account by using the elastic constant matrix whose components are obtained by the tensor transformation from the Cartesian coordinates coincident with the crystallographic axes to those corresponding to a pulling direction [9,101. 2.2. Dislocation density parameters
103
=
o’tOt/l2o-CR5,
2.3. Analysis procedure First of all, the temperature distribution and shape of a silicon bulk single crystal were ob-
Jordan et a!. [121assumed that the dislocation density is proportional to the glide strain, which is caused when the resolved shear stress 0R5 exceeds a certain critical value, the so-called critical resolved shear stress ~CRS’ The slip system in a silicon single crystal is associated with the {111}, (110> slip system which represents twelve permissible glide operations. Taking account of twelve distinct slip systems, Jordan et al. [121 proposed the following dislocation density parameters ~ and o1n1: 12
E I(°~)J~
0ex
(1)
i—i
12
E
0tL)t
i=
(2)
(°Rs)~I’
1
where
a
I
I(0’RS)~
e ~°
~
=
0
UCRS, ftT~~(0R5)l
for
~
> 0CRS’ —
~CRS
b
Fig. 1. Finite element mesh of a single crystal with 6 inch diameter: (a) termination of radial growth steady growth (t = 607 mm).(t
=
104 mm); (b)
104
N. Miyazaki et al.
/ Thermal stress analysis of Si bulk single cn’stal during CZ growth
tamed from a computer program for a transient heat conduction analysis [111, and then a thermal stress analysis was performed, taking account of elastic anisotropy and temperature dependence
a
of material properties [9,10]. Finally, the stress components obtained from the thermal stress analysis were converted into the dislocation density parameters shown in eqs. (1) to (3).
‘‘—‘‘—i’—’~~,’’’’:’~’’r’’’~J’’’’ ~ 8dyn/cm’ ~~ max:i.83~~1O8dyn/Cm2 max:2.g5xlO min:O.O min:O.O A :2.70’lt7dyn/crn2 A :2.7O’~1O7dyn/Cm2
aT”20K
~
rnax:4.85”lO8dyn/cm2 rnin:4.87’lO6dyn/cm2 A :2.70”lO7dyn/crn2
-
max:3.42olD0dyn/cm2 riiln:7.O8olO6dyn/cm2 A :2.7Oo1O~dyn/cm2
400~.
30.0
c 20.0
A
-IV 0
A
0.0
100
A
_____________________ -100 0.0 10.0
cm
cm
Temp.
a~ [001]
—~—‘———i—’——’—’—’——i—’—’——~~1’’’’1’’’’r’’~’
~T—2O — K
ill’ 0.0
100
-I,, .10.0
0.0
A
III.
cm
L:,—,—. ~
~‘‘l’’’’I’’’’I’~’
0.0
0.0
_____________________ 10.0 0.0 10.0
cm
[111]
0~
max:1.8txlO8dyn/crn2 z :2.70x]O7dyn/cm2
30.0
A
~tot
cm
[001]
[111]
0iot
~
nin:O.O max:1.19o1O~dyn/cm2 A :2.7O~’1O7dyn/Cm2
nin:6.38’<1O5dy~/crn2 max:3.tOxlO8dyn/cm2 A :2.7U~1O7dyn/cm~
min:6.96xlQ6dyn/cm2 nax:2.92x10°dyn/cm2 ~ :2.70”lO7dyn/Crfl2
400
30,0~. .
~______ -10.0
0.0
cm Temp.
10.0
-10.0 I~~I
0.0
10.0
T~~T T~I -10.0
cm o~ [001]
0.0
10.0
-10.0
cm [111]
0.0
cm
o~
o~e~[001]
10.0
T~IT
-10.0
0.0
10.0
cm [111]
0t~i
Fig. 2. Contour lines of temperature and the dislocation density parameters o~ and a(a) t
=
104 mm; (b)
t
=
242 mm (c)
t
104 for a single crystal with 6 inch diameter: 607 mm.
N. Miyazaki et al. C
—
51—20K
~
~
500A
:
Thermal stress analysis of Si bulk single crystal during CZ growth
2 min:0.0 max:3.37~<1O8dyfl/crn A :2.70”10~dyn/cm2
max:2.O7xlO0dyn/cm2 min:0.O A :2.70”107dynfcm~
A
A
~
u~.0
/
max:4.79olO8dyn/Cm2 min:L77”lO6dyn/cm2 s :2.70x~07dyn/cm2
105 rnax:3.3lolOOdyn/cm2 min:3.52~106dyn/crn2 A :2.7OolO7dyn/crn2
A
/ N
A
•
~
~
1,1,1 .100
~ 00
I ,~,,,, 1,1,1 100 -10.0
cm Temp.
I. .~.1. ~ 10.0
0.0
.10.0
cm a-~ [001]
I.,II,,,,I 0.0
,,,, 100
III -10.0
II -0.0
cm
Clii
o~ [111J
,,~,,,,, I ~ 400 .10.0
0.0
10.0
cm
o 4~ [001]
0404
[111]
Fig. 2 (continued).
3. Results and discussion The thermal stress analysis of a silicon bulk single crystal with a diameter of 6 or 8 inches during the CZ growth was performed in the cases of the [0011 and [111] pulling directions. Examples of finite element mesh are shown in fig. 1 for two cases, i.e., the termination of radial growth (t la) crystal and steady (t 607 mm 104 ) of mm) a bulk(fig. single with growth 6 inch diameter (fig. Ib). The 20-noded isoparametric three-dimensional solid elements were used in the analysis. The following were the geometrical parameters of a CZ furnace and the operational parameters used in the transient heat conduction analysis to obtain the temperature distribution and shape of a silicon bulk single crystal: inner diameter of a crucible 44 cm, crucible height 34.3 =
=
.
=
=
diameter of 6 inches and 142 kW for a bulk single crystal with a diameter of 8 inches. The heat conduction analysis was performed by fixing heater power and using heat transfer coefficients of rotating bodies estimated from the crucible rotation and the crystal rotation. The material properties required in the therma! stress analysis are2)the elastic constants C11, [14], the thermal expanC12 and C44 (dyn/cm sion coefficient a (K-’) [15], and the critical resolved shear stress 0~CRS (dyn/cm2) [16]. The following were employed in the present analysis: C, 2 1=1.6564x10’ xexp(—9.4 x 105(T— 298.15)), (4a) C 0 6394>< 1012 =
I2
xexp(_9.8 x 105(T— 298.15)),
44 cm, temperature of a furnace wall 300 K, charge of silicon 60 kg,23 crucible —5crystal rpm, crystal rotation rpm for arotation bulk single with a diameter of 6 inches and 16rpm for a bulk single crystal with a diameter of 8 inches, heater power 143 kW for a bulk single crystal with a =
=
=
=
=
(4b)
C a
=
=
0.795 1 x
1012
5(T— 298.15)), (4c) xexp(—8.3 3.725 X i0~ x 10 x [1.0— exp(—5.88 X 103(T— 124.0))] +
5.84
X
10~°T,
(5)
106
N. Miyazaki et al.
a
-
AT2OK
‘
/
Thermal stress analysis of Si bulk single crystal during CZ growth
8dyn/cm2 max:1.02x]O min:0.0 A :2.70”]0~dyn/cm2
max:1.2golt8dyn/cm2 rnin:0.0 A :2.7OolO7dyn/cm2
,—,—,—,—T,,—,—,—, ~
~
oo.o
~
max:3.48x]O0dyn/cm2 mm:] .34o]O7dyn/cm2 A :2.7OolO7dyn/cm2
max:3.92x10°dyn/cni2 min:].47olO7dyn/cm2 A :2.7Oxlo7dynfcm2
400
30.0
20.0
A
___ ‘10.0
I
I
I 00
lU 0
,...!..
I .10.0 ,
0.0
cm Temp.
I
I 0.0
I,
I,,, .100
10.0
cm
AT2OK
I 00
I,,
1.1 0.0
I .100
cm a-~ [111]
[001]
~o~r
I
max:1.94x]O0dyn/cm2 min:0.0 s :2.70x]O7dyn/cm2
. ~,
_ I
00
I,,,., 100
I,,, .10.0
cm ae~~[001]
~ max:1.27x]Q8dyn/cm2 min0.0 A :2.7OxlO7dyrl/cm2
I,
I 0.0
I,
I,. 100
cm ~
max:3.90x]O8dyn/cm2 min:9.28xlO5dyn/cm2 A :2.70~10~dyn/cm2
[111]
max:3.OOolO0dyn/cm’ min:1 .43U]O5dyn/cm2 A :2.70”lO7dyn/cm2
300
~ .10.0
~ 0.0
00.0
~ -100
0,0
cm Temp.
10.0
.100
cm ~
0.0
10.0
cm [001]
0~.r
I,. 10.0
0.0
10.0
~ .10.0
0.0
cm
[111]
at~~[001]
10,0
cm
~
[111]
Fig. 3. Contour lines of temperature and the dislocation density parameters a(a) t
=
104 mm; (b)t
=
242 mm; (c)
t
0~ and a-,,1 for a single crystal with 8 inch diameter: 607 mm.
=
N. Miyazaki et a!.
_
C ‘
40.0
30.0
Thermal stress analysis of Si bulk single crystal during CZ growth
0dyn/cm2 min:0.0 max:2.37olO A :2.7OolO7dyn/cm2
sT=20K 300
/
~
max:].37’~108dyn/Cfl1~
min:0.0 s :2.7OolO7dyn/Cm2
~
max:3.8OolO0dyn/cm2
~ min:4.89o100dyn/croi~ A :2.70o~07dyn/cm2
:
-
107
A,
max:2.93U]Osdyn/cm2
min:6.66olQ6dyn/cmO A :2.70”lO7dyn/cni2
_____________________
E
______
20.0’
I
_
10.0
00
I ,,?,,,1,,1,,l .100 00
.10.0
III -100
,,,,,,
100
II,.’, 10.0
0.0
cm
II -100
cm
Temp.
II,,, 0.0
I,,,,,, 10,0
1,1,1 .10.0
cm
I,
I,,,,, 10.0
0.0
1,1,1 .100
cm
I., 0.0
I,,. 10.0
cm
o 0,0
[001]
a~ [111]
0101
[001]
04~4
[111]
Fig. 3 (continued).
4/T). tTCRs
=
exp(10.55
+
(6)
1.0147 x 10
The unit2.ofThe temperature is K the and 1temperature MPa 10 followingT are dyn/cmof experimental data which were used to ranges obtain eqs. (4)—(6): =
300 to 1000 K,
C 1~and crCRS, have been measured at a temperature much lower than the melting point of silicon, 1685 K. In the present we required material properties nearanalysis, the melting point, so the we used the extrapolated values of eqs. (4)—(6) beyond the above temperature range.
for C 41,
Figs. 2a—2c show the results of a bulk single
300 to 1500 K, for a, 0~CRS~ 300 to 1288 K, for As shown here, the material properties, especially
crystal with a diameter of for three tour lines of temperature T 6andinches two dislocation distinct growth steps. In these figures, the condensity parameters
0~ex
and ~
in a longitudinal
2.0
2.0
b
a 1.0
1.0 E
dh variable
0.0
E ~0.0
constant
~on~region
~ariable
-‘I:,
—1.0
constant region
_
region
—1.0
—2.0
—20 0
2
4
6 (hour)
8
10
0
2
4
6 (hour)
8
10
Fig. 4. Time variation of the shape of the crystal—melt interface: (a) single crystal with 6 inch diameter; (b) single crystal with 8 inch diameter.
108
N. Miyazaki et a!.
/ Thermal stress analysis of Si bulk single crystal during
section are shown for both the [001] and the [111] pulling direction. The same results are shown in figs. 3a—3c for a bulk single crystal with a diameter of 8 inches. The arrows depicted in these figures indicate the locations where the maximum values of °ex and ~ occur. The maximum values occur either at the center of the crystal—melt interface or at the periphery near the crystal—melt interface, Here we discuss the relation between the dislocation density parameters and the shape of the crystal—melt interface, which is characterized by the maximum difference of height between the peripheral and central regions at the bottom of a x10B
bulk single crystal dh. The sign of dh changes from the negative to the positive when the crystal—melt interface changes from a convex shape to a concave one. Figs. 4a and 4b show the dh versus time curves for single crystals with 6 and 8 inch diameters. The dh versus time curve can be divided into a variable region and a constant region, as shown in these figures. In the variable region, a single crystal grows in the radial direction for about the first 2 h, and then it changes its interface shape for the next 4 h, keeping its radius constant. In the constant region, a bulk single crystal continues to grow steadily, keeping both its interface shape and its radius constant. x10B
a
6.0
b
Oex Otat
— ——
‘V
~
_________________________________
6.0
——
CZ growth
1
1
I
~
_____
00 ~Q4~~1~1__11
1,1,,,,I,l~,4,II,II
0
2
4
6 (hour)
8
10
0
8
~
2
4
6 (hour)
~
8
10
xloe
x10 6.0
C
—
I
~
4.0
——
4
0 ex Otot
6.0
d
4.0
—
2.0
2.0
o~o 0
-~‘~—‘‘‘‘‘
2
00
I,44I,I,I
,,~
4
~ ex Otot
——
6 (hour)
8
10
Fig. 5. Time variations of the dislocation density parameters
~
0
2
4
6 (hour)
8
10
a-
00 and a-101: (a) single crystal with 6 inch diameter pulled in the 1001] direction; (b) single crystal with 6 inch diameter pulled in the [111] direction; (c) single crystal with 8 inch diameter pulled in the [001] direction; (d) single crystal with 8 inch diameter inches pulled in the [111]direction.
N. Miyazaki et aL
/
Thermal stress analysis of Si bulk single crystal during CZ growth
The degree of concave shape is larger in a single crystal with 6 inch diameter than in a single crystal with 8 inch diameter, since they have almost the same value of dh, i.e. about 1 cm. Next we show the time variation of the maximum of each dislocation density parameter. The time variations of the maximum values of 0’ex and o~ are depicted in figs. 5a—5d, and those of °‘niex are shown in figs. 6a and 6b. In these figures, the results of single crystals with 6 and 8 inch diameter are shown in the cases of the [001] and [111] pulling directions. Comparing these figures with figs. 4a and 4b we can find the relation between the time variations of the dislocation density parameters and that of the shape of the crystal—melt interface. The dislocation density parameters vary greatly during the radial growth, and they show the minimum when the shape of the crystal—melt interface becomes fiat. Speaking roughly, the dislocation density parameters are proportional to the absolute value of dh. This is due to the following reason. The radial gradient of the temperature in a single crystal near the interface decreases when the shape of the interface becomes flat, because the shape of the interface represents the contour line of the melting point, The thermal stress and the dislocation density parameters decrease as a result of the decrease in the radial gradient of the temperature. After terminating radial growth, the dislocation density parameters increase with the increase of concav-
109
ity of the interface, and they show nearly constant values in the constant region shown in figs. 4a and 4b. In the constant region, a single crystal with 8 inch diameter shows smaller dislocation density parameters than that with 6 inch diameter, because both crystals have almost the same value of dh in the constant region and then the concavity of the former is shallower than that of the latter. A single crystal with 8 inch diameter has lower dislocation density parameters than that with 6 inch diameter, as shown in figs. 5 and 6. This may be due to the fact that a single crystal with 8 inch diameter has smaller temperature gradient because of larger radiation from a crucible wall which results from smaller distance between a crystal surface and a crucible wall than a single crystal with 6 inch diameter. Comparing the results of the [0011and [111] pulling directions, the result of the [111] pulling direction is about 60% to 70% of that of the [001] pulling direction. The time variations of the location where the maximum of 0ex or ~~ 00arises are shown in figs. 7a—7d by giving the radius and height of its location. The maximum values of ~ and O~1oi arise at the periphery near the crystal—melt interface during the radial growth. After terminating the radial growth, they arise either at the periphery near the interface or at the center near the interface, depending on the pulling direction and the size of a single crystal.
11
4
TJ4~ 111411T11T1111111T14111411!111111
11T1
4
TTTjIiTfuui1uuijlri1iiIjlhljllIjIIl1luI~lli
—
———
0
Ciii]
—cool]
b
—cool] —
———
Clii)
-
0111111
(hour)
(hour)
Fig. 6. Time variation of the dislocation density parameter a-moo: (a) single crystal with 6 inch diameter; (b) single crystal with 8 inch diameter.
111)
N. Miyazaki et al.
3o.0
~
/ Thermal stress analysis of Si bulk single crystal during CZ growth
a
30.0
‘‘iii’’’
,,,II,,lii,i,I,,
..
0 cx 0 ex
—
20.0
radius height
b
Otot radius —
—
20.0
0 tot height
E
~10.0
-
—to.c
-
E ~-10.0
~
0.0
0 ex 0 ex
-
——— — —
2
4
6
___________________________________________________
8
10
0
2
4
—
200
——— — —
Ocx 0 cx
radius height
Otot 0 tot height radius
~“i’i
~‘::
8
10
2
4
6 (hour)
d
Ocx 0 cx
20.0
——— — —
radius height
0Otot tot height radius,
~1::
l1lIli4illiIil~IlillllillllI1lilllll~lllli~
~
1.1_I_L_Luuluuluulhuulhuhuuuluu,llhhuhhhu
o
6 (hour)
30.0
~
c
—10.0
Qtot radius, 0 tot height
-
(hour)
30.0
radius height
0.0
iiili,iIiiiliilI4liI4ilIliillilIl~iliilIlll
0
~
8
10
0
2
4
6 (hour)
8
10
Fig. 7. Time variations of the locations where the maximum values of the dislocation density parameters a-~
0 and a-,4,1 arise: (a) single crystal with 6 inch diameter pulled in the [001] direction; (b) single crystal with 6 inch diameter pulled in the [111] direction: (c) single crystal with 8 inch diameter pulled in the [001] direction; (d) single crystal with 8 inch diameter pulled in the [Ill] direction.
4. Concluding remarks The thermal stress analysis of a silicon bulk single crystal with 6 or 8 inch diameter was performed and the relation between the dislocation
density parameters and the shape of the crystal— melt interface was discussed. The correlation between the values of the dislocation density parameters and the shape of the interface was confirmed by the present analysis.
N. Miyazaki et al.
/
Thermal stress analysis of Si bulk single crystal during CZ growth
Dislocations are expected to be formed in the present study, because the dislocation density parameters o~ and a-,~ are, respectively, larger than zero and unity. On the other hand, it is known that dislocation-free bulk single crystals of silicon can usually be produced by the CZ growth in the silicon industry. One of the reasons for this paradox may he the uncertainty of the material properties of silicon near the melting point. In the present analysis, extrapolated values were used as the material properties near the melting point. We are not sure if these values are reliable. Experimental or theoretical studies should be done to determine the material properties of silicon near the melting point. One of the most promising ways may be molecular dynamics. In order to make a quantitative comparison between analyses and experiments, we should perform the analysis using reliable material properties near the melting point.
Acknowledgement This study was partially supported by a Grantin-Aid for Scientific Research from the Ministry of Education, Science and Culture.
Ill
References [1] AS. Jordan, R. Caruso and AR. Von Neida, Bell System Tech. J. 59 (1980) 593. [2] M.J. Duseaux, J. Crystal Growth 61(1983) 576. [3] N. Kobayashi and T. twaki, J. Crystal Growth 73 (1985) 96. [4]G. Szabo, J. Crystal Growth 73 (1985) 131. [51J.C. Lambropoulos and C.N. Delametter, J Crystal Growth 92 (1988) 390. [6]CE. Schvezov. IV. Samarasekera and F. Weinberg. J. Crystal Growth 92 (1988) 479. [7] S. Motakef, J. Crystal Growth 96 (1989) 201. [8]F. Dupret, P. Nicodeme and Y Ryckmans, J. Crystal Growth 97 (1989) 162. [91N. Miyazaki, S. Hagihara and T. Munakata, J. Crystal Growth 106 (1990) 149. [101N. Miyazaki, H. Uchida, S. Hagihara and T. Munakata, J. Crystal Growth 113 (1991) 227. [11]W. Nakayama and H. Masaki, in: Heat Transfer — 1986, Vol 4 (Hemisphere, 1986) p. 1755 [121 AS. Jordan, R. Caruso, AR. Von Neida and J.W. Nielsen, J. AppI. Phys. 52(1981) 3331. [13] H. Ohya, Proc. Japan Soc. Mech. Eng. No. 910-17, Vol. A (1991) 312. 114] INSPEC. Properties of Silicon, EMIS Datareviews Series, No. 4 (1988) 14 [15] Y. Okada andY. Tokumaru, J. AppI. Phys. 56(1984)314. [16]AS. Jordan, R. Caruso and AR. Von Neida, J. Crystal Growth 79 (1986) 243.
j
~
CRYSTAL
GROW T H
Thermal stress analysis of silicon bulk single crystal during Czochralski growth N. Miyazaki, H. Uchida, T. Munakata Department of Chemical Engineering, Faculty of Engineering, Kyushu UniLersity. 6-10-1 Hakozak,, Higashi-ku, Fukuoka-shi, Fukuoka-ken 812, Japan
K. Fujioka Mechanical Engineering Research Laboratory, Hitachi, Ltd. 502 Kandatsu-machi, Tsuchiura-shi, Ibaraki-ken 300, Japan
and Y. Sugino Kofu Works, Hitachi, Ltd., Nishiyahata, Ryuoh-cho, Nakakoma-gun, Yamanashi-ken 400-01, Japan Received 15 April 1992; manuscript received in final form 21 June 1992
The thermal stress analysis of a silicon bulk single crystal with a diameter of 6 or 8 inches is performed in the cases of the [001] and [1111pulling directions by using a three-dimensional finite element program developed for calculating thermal stress in a bulk single crystal during the Czochralski growth. Elastic anisotropy and temperature dependence of material properties are taken into account in this program. The temperature distribution and shape of a silicon bulk single crystal which are required for the thermal stress analysis are obtained from a computer program for a transient heat conduction analysis which is specialized for the Czochralski growth. The stress components obtained from the thermal stress analysis are converted into the parameters related with dislocation density. The time variations of these parameters are shown in this paper. The relation between these parameters and the shape of the crystal—melt interface is discussed.
1. Introduction The dislocations formed in bulk single crystals during the Czochralski (CZ) growth have adverse effects on the performance of electronic devices, It is known that the thermal stress during the CZ growth plays an important role in the formation of dislocations, Although a lot of studies [1—8] have been performed concerning the thermal stress analysis of bulk single crystals during the CZ growth, there have been few studies which deal with the thermal stress analysis of such large-sized bulk single crystals as industrial use during a long period of the CZ growth. 0022-0248/92/$05.00 © 1992
—
Miyazaki et a!. [9,10] have developed a threedimensional finite element program for the thermal stress analysis of single crystals, which considers the elastic anisotropy and the temperature dependence of the material properties. On the other hand, Nakayama and Masaki [11] have developed a computer program for a transient heat conduction analysis which can simulate the CZ growth process of a bulk single crystal by producing finite elements automatically and by considering heat radiation between the bulk single crystal and its neighboring environment, i.e., melt, heater and radiation shield. In the present study, the thermal stress analysis of a silicon bulk single crystal with a diameter
Elsevier Science Publishers By. All rights reserved
N. Miyazaki et al.
/ Thermal stress analysis
of 6 or 8 inches during a long period of the CZ growth is performed by using both the abovementioned computer programs for a thermal stress analysis and a transient heat conduction analysis. The stress components obtained from the thermal stress analysis are converted into the parameters related with dislocation density, i.e., dislocation density parameters. The time variations of these parameters are shown and the relation between these parameters and the shape of the crystal—melt interface is discussed.
of Si bulk single crystal during CZ growth
The resolved shear stress can be calculated by the transformation of stresses from the Cartesian coordinates used in a thermal stress analysis to a slip system of a single crystal. When 0~ex is equal to zero, no dislocation is expected to be generated. Speaking exactly, 0~tot is not a dislocation density parameter but a measure for a thermal stress effective to glide strain. In addition to o~ and ~ we also used another dislocation density parameter °‘mex [13], which is defined as (3) When °~mexis less than unity, no dislocation is expected to be generated. The dislocation density parameter 0~ex indicates whether or not the thermal stress level in a bulk single crystal exceeds 0~cRs’but does not indicate how many times the thermal stress level is as large as ~CRS’ Therefore, we used another dislocation density parameter 0’mex as a quantitative measure indicating the difference between the thermal stress level and ~
°mex
2. Method of analysis 2.1. Elastic anisotropy In the present analysis, the elastic anisotropy is taken into account by using the elastic constant matrix whose components are obtained by the tensor transformation from the Cartesian coordinates coincident with the crystallographic axes to those corresponding to a pulling direction [9,101. 2.2. Dislocation density parameters
103
=
o’tOt/l2o-CR5,
2.3. Analysis procedure First of all, the temperature distribution and shape of a silicon bulk single crystal were ob-
Jordan et a!. [121assumed that the dislocation density is proportional to the glide strain, which is caused when the resolved shear stress 0R5 exceeds a certain critical value, the so-called critical resolved shear stress ~CRS’ The slip system in a silicon single crystal is associated with the {111}, (110> slip system which represents twelve permissible glide operations. Taking account of twelve distinct slip systems, Jordan et al. [121 proposed the following dislocation density parameters ~ and o1n1: 12
E I(°~)J~
0ex
(1)
i—i
12
E
0tL)t
i=
(2)
(°Rs)~I’
1
where
a
I
I(0’RS)~
e ~°
~
=
0
UCRS, ftT~~(0R5)l
for
~
> 0CRS’ —
~CRS
b
Fig. 1. Finite element mesh of a single crystal with 6 inch diameter: (a) termination of radial growth steady growth (t = 607 mm).(t
=
104 mm); (b)
104
N. Miyazaki et al.
/ Thermal stress analysis of Si bulk single cn’stal during CZ growth
tamed from a computer program for a transient heat conduction analysis [111, and then a thermal stress analysis was performed, taking account of elastic anisotropy and temperature dependence
a
of material properties [9,10]. Finally, the stress components obtained from the thermal stress analysis were converted into the dislocation density parameters shown in eqs. (1) to (3).
‘‘—‘‘—i’—’~~,’’’’:’~’’r’’’~J’’’’ ~ 8dyn/cm’ ~~ max:i.83~~1O8dyn/Cm2 max:2.g5xlO min:O.O min:O.O A :2.70’lt7dyn/crn2 A :2.7O’~1O7dyn/Cm2
aT”20K
~
rnax:4.85”lO8dyn/cm2 rnin:4.87’lO6dyn/cm2 A :2.70”lO7dyn/crn2
-
max:3.42olD0dyn/cm2 riiln:7.O8olO6dyn/cm2 A :2.7Oo1O~dyn/cm2
400~.
30.0
c 20.0
A
-IV 0
A
0.0
100
A
_____________________ -100 0.0 10.0
cm
cm
Temp.
a~ [001]
—~—‘———i—’——’—’—’——i—’—’——~~1’’’’1’’’’r’’~’
~T—2O — K
ill’ 0.0
100
-I,, .10.0
0.0
A
III.
cm
L:,—,—. ~
~‘‘l’’’’I’’’’I’~’
0.0
0.0
_____________________ 10.0 0.0 10.0
cm
[111]
0~
max:1.8txlO8dyn/crn2 z :2.70x]O7dyn/cm2
30.0
A
~tot
cm
[001]
[111]
0iot
~
nin:O.O max:1.19o1O~dyn/cm2 A :2.7O~’1O7dyn/Cm2
nin:6.38’<1O5dy~/crn2 max:3.tOxlO8dyn/cm2 A :2.7U~1O7dyn/cm~
min:6.96xlQ6dyn/cm2 nax:2.92x10°dyn/cm2 ~ :2.70”lO7dyn/Crfl2
400
30,0~. .
~______ -10.0
0.0
cm Temp.
10.0
-10.0 I~~I
0.0
10.0
T~~T T~I -10.0
cm o~ [001]
0.0
10.0
-10.0
cm [111]
0.0
cm
o~
o~e~[001]
10.0
T~IT
-10.0
0.0
10.0
cm [111]
0t~i
Fig. 2. Contour lines of temperature and the dislocation density parameters o~ and a(a) t
=
104 mm; (b)
t
=
242 mm (c)
t
104 for a single crystal with 6 inch diameter: 607 mm.
N. Miyazaki et al. C
—
51—20K
~
~
500A
:
Thermal stress analysis of Si bulk single crystal during CZ growth
2 min:0.0 max:3.37~<1O8dyfl/crn A :2.70”10~dyn/cm2
max:2.O7xlO0dyn/cm2 min:0.O A :2.70”107dynfcm~
A
A
~
u~.0
/
max:4.79olO8dyn/Cm2 min:L77”lO6dyn/cm2 s :2.70x~07dyn/cm2
105 rnax:3.3lolOOdyn/cm2 min:3.52~106dyn/crn2 A :2.7OolO7dyn/crn2
A
/ N
A
•
~
~
1,1,1 .100
~ 00
I ,~,,,, 1,1,1 100 -10.0
cm Temp.
I. .~.1. ~ 10.0
0.0
.10.0
cm a-~ [001]
I.,II,,,,I 0.0
,,,, 100
III -10.0
II -0.0
cm
Clii
o~ [111J
,,~,,,,, I ~ 400 .10.0
0.0
10.0
cm
o 4~ [001]
0404
[111]
Fig. 2 (continued).
3. Results and discussion The thermal stress analysis of a silicon bulk single crystal with a diameter of 6 or 8 inches during the CZ growth was performed in the cases of the [0011 and [111] pulling directions. Examples of finite element mesh are shown in fig. 1 for two cases, i.e., the termination of radial growth (t la) crystal and steady (t 607 mm 104 ) of mm) a bulk(fig. single with growth 6 inch diameter (fig. Ib). The 20-noded isoparametric three-dimensional solid elements were used in the analysis. The following were the geometrical parameters of a CZ furnace and the operational parameters used in the transient heat conduction analysis to obtain the temperature distribution and shape of a silicon bulk single crystal: inner diameter of a crucible 44 cm, crucible height 34.3 =
=
.
=
=
diameter of 6 inches and 142 kW for a bulk single crystal with a diameter of 8 inches. The heat conduction analysis was performed by fixing heater power and using heat transfer coefficients of rotating bodies estimated from the crucible rotation and the crystal rotation. The material properties required in the therma! stress analysis are2)the elastic constants C11, [14], the thermal expanC12 and C44 (dyn/cm sion coefficient a (K-’) [15], and the critical resolved shear stress 0~CRS (dyn/cm2) [16]. The following were employed in the present analysis: C, 2 1=1.6564x10’ xexp(—9.4 x 105(T— 298.15)), (4a) C 0 6394>< 1012 =
I2
xexp(_9.8 x 105(T— 298.15)),
44 cm, temperature of a furnace wall 300 K, charge of silicon 60 kg,23 crucible —5crystal rpm, crystal rotation rpm for arotation bulk single with a diameter of 6 inches and 16rpm for a bulk single crystal with a diameter of 8 inches, heater power 143 kW for a bulk single crystal with a =
=
=
=
=
(4b)
C a
=
=
0.795 1 x
1012
5(T— 298.15)), (4c) xexp(—8.3 3.725 X i0~ x 10 x [1.0— exp(—5.88 X 103(T— 124.0))] +
5.84
X
10~°T,
(5)
106
N. Miyazaki et al.
a
-
AT2OK
‘
/
Thermal stress analysis of Si bulk single crystal during CZ growth
8dyn/cm2 max:1.02x]O min:0.0 A :2.70”]0~dyn/cm2
max:1.2golt8dyn/cm2 rnin:0.0 A :2.7OolO7dyn/cm2
,—,—,—,—T,,—,—,—, ~
~
oo.o
~
max:3.48x]O0dyn/cm2 mm:] .34o]O7dyn/cm2 A :2.7OolO7dyn/cm2
max:3.92x10°dyn/cni2 min:].47olO7dyn/cm2 A :2.7Oxlo7dynfcm2
400
30.0
20.0
A
___ ‘10.0
I
I
I 00
lU 0
,...!..
I .10.0 ,
0.0
cm Temp.
I
I 0.0
I,
I,,, .100
10.0
cm
AT2OK
I 00
I,,
1.1 0.0
I .100
cm a-~ [111]
[001]
~o~r
I
max:1.94x]O0dyn/cm2 min:0.0 s :2.70x]O7dyn/cm2
. ~,
_ I
00
I,,,., 100
I,,, .10.0
cm ae~~[001]
~ max:1.27x]Q8dyn/cm2 min0.0 A :2.7OxlO7dyrl/cm2
I,
I 0.0
I,
I,. 100
cm ~
max:3.90x]O8dyn/cm2 min:9.28xlO5dyn/cm2 A :2.70~10~dyn/cm2
[111]
max:3.OOolO0dyn/cm’ min:1 .43U]O5dyn/cm2 A :2.70”lO7dyn/cm2
300
~ .10.0
~ 0.0
00.0
~ -100
0,0
cm Temp.
10.0
.100
cm ~
0.0
10.0
cm [001]
0~.r
I,. 10.0
0.0
10.0
~ .10.0
0.0
cm
[111]
at~~[001]
10,0
cm
~
[111]
Fig. 3. Contour lines of temperature and the dislocation density parameters a(a) t
=
104 mm; (b)t
=
242 mm; (c)
t
0~ and a-,,1 for a single crystal with 8 inch diameter: 607 mm.
=
N. Miyazaki et a!.
_
C ‘
40.0
30.0
Thermal stress analysis of Si bulk single crystal during CZ growth
0dyn/cm2 min:0.0 max:2.37olO A :2.7OolO7dyn/cm2
sT=20K 300
/
~
max:].37’~108dyn/Cfl1~
min:0.0 s :2.7OolO7dyn/Cm2
~
max:3.8OolO0dyn/cm2
~ min:4.89o100dyn/croi~ A :2.70o~07dyn/cm2
:
-
107
A,
max:2.93U]Osdyn/cm2
min:6.66olQ6dyn/cmO A :2.70”lO7dyn/cni2
_____________________
E
______
20.0’
I
_
10.0
00
I ,,?,,,1,,1,,l .100 00
.10.0
III -100
,,,,,,
100
II,.’, 10.0
0.0
cm
II -100
cm
Temp.
II,,, 0.0
I,,,,,, 10,0
1,1,1 .10.0
cm
I,
I,,,,, 10.0
0.0
1,1,1 .100
cm
I., 0.0
I,,. 10.0
cm
o 0,0
[001]
a~ [111]
0101
[001]
04~4
[111]
Fig. 3 (continued).
4/T). tTCRs
=
exp(10.55
+
(6)
1.0147 x 10
The unit2.ofThe temperature is K the and 1temperature MPa 10 followingT are dyn/cmof experimental data which were used to ranges obtain eqs. (4)—(6): =
300 to 1000 K,
C 1~and crCRS, have been measured at a temperature much lower than the melting point of silicon, 1685 K. In the present we required material properties nearanalysis, the melting point, so the we used the extrapolated values of eqs. (4)—(6) beyond the above temperature range.
for C 41,
Figs. 2a—2c show the results of a bulk single
300 to 1500 K, for a, 0~CRS~ 300 to 1288 K, for As shown here, the material properties, especially
crystal with a diameter of for three tour lines of temperature T 6andinches two dislocation distinct growth steps. In these figures, the condensity parameters
0~ex
and ~
in a longitudinal
2.0
2.0
b
a 1.0
1.0 E
dh variable
0.0
E ~0.0
constant
~on~region
~ariable
-‘I:,
—1.0
constant region
_
region
—1.0
—2.0
—20 0
2
4
6 (hour)
8
10
0
2
4
6 (hour)
8
10
Fig. 4. Time variation of the shape of the crystal—melt interface: (a) single crystal with 6 inch diameter; (b) single crystal with 8 inch diameter.
108
N. Miyazaki et a!.
/ Thermal stress analysis of Si bulk single crystal during
section are shown for both the [001] and the [111] pulling direction. The same results are shown in figs. 3a—3c for a bulk single crystal with a diameter of 8 inches. The arrows depicted in these figures indicate the locations where the maximum values of °ex and ~ occur. The maximum values occur either at the center of the crystal—melt interface or at the periphery near the crystal—melt interface, Here we discuss the relation between the dislocation density parameters and the shape of the crystal—melt interface, which is characterized by the maximum difference of height between the peripheral and central regions at the bottom of a x10B
bulk single crystal dh. The sign of dh changes from the negative to the positive when the crystal—melt interface changes from a convex shape to a concave one. Figs. 4a and 4b show the dh versus time curves for single crystals with 6 and 8 inch diameters. The dh versus time curve can be divided into a variable region and a constant region, as shown in these figures. In the variable region, a single crystal grows in the radial direction for about the first 2 h, and then it changes its interface shape for the next 4 h, keeping its radius constant. In the constant region, a bulk single crystal continues to grow steadily, keeping both its interface shape and its radius constant. x10B
a
6.0
b
Oex Otat
— ——
‘V
~
_________________________________
6.0
——
CZ growth
1
1
I
~
_____
00 ~Q4~~1~1__11
1,1,,,,I,l~,4,II,II
0
2
4
6 (hour)
8
10
0
8
~
2
4
6 (hour)
~
8
10
xloe
x10 6.0
C
—
I
~
4.0
——
4
0 ex Otot
6.0
d
4.0
—
2.0
2.0
o~o 0
-~‘~—‘‘‘‘‘
2
00
I,44I,I,I
,,~
4
~ ex Otot
——
6 (hour)
8
10
Fig. 5. Time variations of the dislocation density parameters
~
0
2
4
6 (hour)
8
10
a-
00 and a-101: (a) single crystal with 6 inch diameter pulled in the 1001] direction; (b) single crystal with 6 inch diameter pulled in the [111] direction; (c) single crystal with 8 inch diameter pulled in the [001] direction; (d) single crystal with 8 inch diameter inches pulled in the [111]direction.
N. Miyazaki et aL
/
Thermal stress analysis of Si bulk single crystal during CZ growth
The degree of concave shape is larger in a single crystal with 6 inch diameter than in a single crystal with 8 inch diameter, since they have almost the same value of dh, i.e. about 1 cm. Next we show the time variation of the maximum of each dislocation density parameter. The time variations of the maximum values of 0’ex and o~ are depicted in figs. 5a—5d, and those of °‘niex are shown in figs. 6a and 6b. In these figures, the results of single crystals with 6 and 8 inch diameter are shown in the cases of the [001] and [111] pulling directions. Comparing these figures with figs. 4a and 4b we can find the relation between the time variations of the dislocation density parameters and that of the shape of the crystal—melt interface. The dislocation density parameters vary greatly during the radial growth, and they show the minimum when the shape of the crystal—melt interface becomes fiat. Speaking roughly, the dislocation density parameters are proportional to the absolute value of dh. This is due to the following reason. The radial gradient of the temperature in a single crystal near the interface decreases when the shape of the interface becomes flat, because the shape of the interface represents the contour line of the melting point, The thermal stress and the dislocation density parameters decrease as a result of the decrease in the radial gradient of the temperature. After terminating radial growth, the dislocation density parameters increase with the increase of concav-
109
ity of the interface, and they show nearly constant values in the constant region shown in figs. 4a and 4b. In the constant region, a single crystal with 8 inch diameter shows smaller dislocation density parameters than that with 6 inch diameter, because both crystals have almost the same value of dh in the constant region and then the concavity of the former is shallower than that of the latter. A single crystal with 8 inch diameter has lower dislocation density parameters than that with 6 inch diameter, as shown in figs. 5 and 6. This may be due to the fact that a single crystal with 8 inch diameter has smaller temperature gradient because of larger radiation from a crucible wall which results from smaller distance between a crystal surface and a crucible wall than a single crystal with 6 inch diameter. Comparing the results of the [0011and [111] pulling directions, the result of the [111] pulling direction is about 60% to 70% of that of the [001] pulling direction. The time variations of the location where the maximum of 0ex or ~~ 00arises are shown in figs. 7a—7d by giving the radius and height of its location. The maximum values of ~ and O~1oi arise at the periphery near the crystal—melt interface during the radial growth. After terminating the radial growth, they arise either at the periphery near the interface or at the center near the interface, depending on the pulling direction and the size of a single crystal.
11
4
TJ4~ 111411T11T1111111T14111411!111111
11T1
4
TTTjIiTfuui1uuijlri1iiIjlhljllIjIIl1luI~lli
—
———
0
Ciii]
—cool]
b
—cool] —
———
Clii)
-
0111111
(hour)
(hour)
Fig. 6. Time variation of the dislocation density parameter a-moo: (a) single crystal with 6 inch diameter; (b) single crystal with 8 inch diameter.
111)
N. Miyazaki et al.
3o.0
~
/ Thermal stress analysis of Si bulk single crystal during CZ growth
a
30.0
‘‘iii’’’
,,,II,,lii,i,I,,
..
0 cx 0 ex
—
20.0
radius height
b
Otot radius —
—
20.0
0 tot height
E
~10.0
-
—to.c
-
E ~-10.0
~
0.0
0 ex 0 ex
-
——— — —
2
4
6
___________________________________________________
8
10
0
2
4
—
200
——— — —
Ocx 0 cx
radius height
Otot 0 tot height radius
~“i’i
~‘::
8
10
2
4
6 (hour)
d
Ocx 0 cx
20.0
——— — —
radius height
0Otot tot height radius,
~1::
l1lIli4illiIil~IlillllillllI1lilllll~lllli~
~
1.1_I_L_Luuluuluulhuulhuhuuuluu,llhhuhhhu
o
6 (hour)
30.0
~
c
—10.0
Qtot radius, 0 tot height
-
(hour)
30.0
radius height
0.0
iiili,iIiiiliilI4liI4ilIliillilIl~iliilIlll
0
~
8
10
0
2
4
6 (hour)
8
10
Fig. 7. Time variations of the locations where the maximum values of the dislocation density parameters a-~
0 and a-,4,1 arise: (a) single crystal with 6 inch diameter pulled in the [001] direction; (b) single crystal with 6 inch diameter pulled in the [111] direction: (c) single crystal with 8 inch diameter pulled in the [001] direction; (d) single crystal with 8 inch diameter pulled in the [Ill] direction.
4. Concluding remarks The thermal stress analysis of a silicon bulk single crystal with 6 or 8 inch diameter was performed and the relation between the dislocation
density parameters and the shape of the crystal— melt interface was discussed. The correlation between the values of the dislocation density parameters and the shape of the interface was confirmed by the present analysis.
N. Miyazaki et al.
/
Thermal stress analysis of Si bulk single crystal during CZ growth
Dislocations are expected to be formed in the present study, because the dislocation density parameters o~ and a-,~ are, respectively, larger than zero and unity. On the other hand, it is known that dislocation-free bulk single crystals of silicon can usually be produced by the CZ growth in the silicon industry. One of the reasons for this paradox may he the uncertainty of the material properties of silicon near the melting point. In the present analysis, extrapolated values were used as the material properties near the melting point. We are not sure if these values are reliable. Experimental or theoretical studies should be done to determine the material properties of silicon near the melting point. One of the most promising ways may be molecular dynamics. In order to make a quantitative comparison between analyses and experiments, we should perform the analysis using reliable material properties near the melting point.
Acknowledgement This study was partially supported by a Grantin-Aid for Scientific Research from the Ministry of Education, Science and Culture.
Ill
References [1] AS. Jordan, R. Caruso and AR. Von Neida, Bell System Tech. J. 59 (1980) 593. [2] M.J. Duseaux, J. Crystal Growth 61(1983) 576. [3] N. Kobayashi and T. twaki, J. Crystal Growth 73 (1985) 96. [4]G. Szabo, J. Crystal Growth 73 (1985) 131. [51J.C. Lambropoulos and C.N. Delametter, J Crystal Growth 92 (1988) 390. [6]CE. Schvezov. IV. Samarasekera and F. Weinberg. J. Crystal Growth 92 (1988) 479. [7] S. Motakef, J. Crystal Growth 96 (1989) 201. [8]F. Dupret, P. Nicodeme and Y Ryckmans, J. Crystal Growth 97 (1989) 162. [91N. Miyazaki, S. Hagihara and T. Munakata, J. Crystal Growth 106 (1990) 149. [101N. Miyazaki, H. Uchida, S. Hagihara and T. Munakata, J. Crystal Growth 113 (1991) 227. [11]W. Nakayama and H. Masaki, in: Heat Transfer — 1986, Vol 4 (Hemisphere, 1986) p. 1755 [121 AS. Jordan, R. Caruso, AR. Von Neida and J.W. Nielsen, J. AppI. Phys. 52(1981) 3331. [13] H. Ohya, Proc. Japan Soc. Mech. Eng. No. 910-17, Vol. A (1991) 312. 114] INSPEC. Properties of Silicon, EMIS Datareviews Series, No. 4 (1988) 14 [15] Y. Okada andY. Tokumaru, J. AppI. Phys. 56(1984)314. [16]AS. Jordan, R. Caruso and AR. Von Neida, J. Crystal Growth 79 (1986) 243.