Modelling the dynamics and control design for Czochralski, Liquid Encapsulated Czochralski and Floating Zone processes

Available online at www.sciencedirect.com

Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121 www.elsevier.com/locate/pcrysgrow

Modelling the dynamics and control design for Czochralski, Liquid Encapsulated Czochralski and Floating Zone processes Gennadii Satunkin a,b,* a

ITV Inc., Intellect.Technologies.VIDEO, Moscow-Chernogolovka, Russia b AxxonSoft Inc, Moscow-Chernogolovka, Russia

Abstract The need for high quality crystals constantly rise. It is especially obvious in connection with the evolution of electronics and optoelectronics. At the present time all the basic methods for crystal growing are known. So the question arises e what further developments are needed to create further advances? Without doubt first of all it is necessary to speak about perfecting crystal growing equipment. Perfecting modern equipment will enhance and accommodate the results of our understanding of the crystallization physics and provide solutions to the various physical tasks at the atomic and macroscopic levels. Each new step in the process of perfecting the technology and production processes demands large intellectual and material inputs. The continual updating of pullers requires constructive solutions and control systems. Mathematical modelling of the methods of crystallization enables one to more rapidly create the software for the digital systems which are a feature of the latest achievements of physics, that is IT engineering and the modern theory of automatic control. Here we consider the problem of mathematical modelling of crystallization processes from the melt by Czochralski, Liquid Encapsulated Czochralski and Floating Zone methods based on linearization of three conservation laws: the heat, mass and the growth angle constancy is reviewed in depth. Special attention is given to the problem of the dynamical analysis of these processes in open and closed states and to the synthesis of the digital control of crystal diameter for the weight technique. The main problems discussed involve the determination of the parameters required for calculating such control systems together with the use of the multichannel parametric PID regulator involving the state variable observer concept. In addition

* Corresponding author. 142 432 Chernogolovka, Centralnaja 18-55, Moscow District, Russia. Tel.: þ7 496 52 2 33 73, þ79104798649 (mobile). E-mail address: [email protected]. URL: http://www.itv.ru, http://www.axxonsoft.com. 0960-8974/$ - see front matter
2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.pcrysgrow.2010.05.001

2

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

briefly considered are the problems of digital filtering of the measurement noise based on multidimentional Kalman filters and the determination of mechanical stability limits for static menisci.
2010 Elsevier Ltd. All rights reserved. Keywords: A1. Mathematical modelling; A1. Modern control theory; A2. Automated melt growth; A2. Czochralski Growth; A2. Liquid Encapsulation Czochralski (LEC); A2. Floating zone technique; B1. Semiconductors; B1. oxides

1. Introduction. The crystal surface formation in the Czochralski method The problem of lateral crystal surface formation in the vicinity of three-phase line (TPL) consists in the determination of the orientation of a new external solid surface element for the melt surface in a defined position. In the Cz method the meniscus surface orientation depends on the shape and the height h of the TPL L (Fig. 1) and in general is correctly determined by the relationship dr(t)/dl(t)¼tg4s(t) or dl(t) ¼ Vcr(t)dt, consequently: dr_ ðtÞ ¼ Vcr ðtÞtg4s ðtÞ

ð1:1Þ

The modern Voronkov theory [1e7] of surface formation for growing crystal takes into account the crystalline surface anisotropy, the surface mass transfer and is based on four classical results: 1. the Wulff theorem about the equilibrium crystal shape [8], 2. the Herring theorems about free surface energies and equilibrium of multiple crystalline anisotropic surfaces at a junction line [9,10], 3. the Hoffman-Cahn geometrical construction for the equilibrium state of crystalline surfaces at a junction line [11], 4. and also the Bolling and Tiller analysis of equilibrium shape of a solid phase during crystallization from the melt [12].

Fig. 1. The Crystallization rate Vcr(t) differs from the crystal pulling rate V0 since the crystallization front displacement is dependent on the thermal conditions in the growth zone.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

3

Let’s very briefly consider these positions. (1) The equilibrium crystal form is determined RR by the minimum of the first variation of the crystal free surface energy functional d sðnÞds ¼ min (with the restriction, that the volume of solid phase is constant). The task involves the solution that represents the envelope family of planes P(n) which satisfy to the equation (r$n) ¼ s(n), where r is a vector from center O to the point on the plane P(n) [4,5,8]. The geometric solution of this equation is based on using the Wulff approach and the polar diagram. The equilibrium crystal surface coincides with the envelope G (Wulff figure) of planes Pi(ni) (Fig. 2). In the Cz method it is possible, that the two crystalline interfaces (front and lateral surface) may be singular, vicinal or weakly anisotropic. These surfaces are classified as stable, metastable and labile [4,5,8,13]. For RR the stable surfaces the second variation of the free surface energy functional is positive d2 sðnÞds > 0. The stable orientation n must have a positive radius of curvature s(n) þ s"(n) > 0 at the corresponding point on the Wulff figure. Here s"(n) e is the second angular differential of the free energy s(n) [8]. The singular orientations correspond to planar parts on the envelope G (Fig. 2) since the angular differential of the polar diagram in this case is s"(n) ¼ 0. The length of such planar parts is determined by the magnitude of the jump s’(n0) [4,5,13]. In Fig. 2 the metastable orientations are represented by the dotted lines. The labile orientations, with concave (outside) portions of G, are not shown. Let’s use below for the vectors of the free surface energies sLG, sSG and sSL notations sm, sc and sf and also define ni and ei as the normal and tangent vectors to the surfaces i ¼ m, c, f. (2) According to the Herring theorem, the general equilibrium condition for i intersecting surfaces (i3) is determined by the variation of a free surface a small deflection Penergy for = of the anisotropic surfaces and is equal to zero dF ¼  ðsi e þ si nÞds ¼ 0 [9]. It is possible to rewrite this condition like the vector balance: X ðsi e þ s0i nÞ ¼ 0 ð1:2Þ i

(3) The vector thermodynamic approach of Hoffman-Cahn [11], based on a graphical solution of the Herring equation set (1.2), allows one to investigate the problem of the crystal

Fig. 2. The envelope G of planes Pi(ni) is the section of the spatial Wulff figure. The stable crystalline surfaces Pi correspond to the convex (outside) portion of the envelope. Dotted lines show the metastable orientations (surfaces) on the Wulff figure G. Planar parts on G correspond to the singular orientations (minimum on the polar diagram F).

4

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

surface formation in the general case, i.e. the presence of singular facets at the TPL [2e7]. The graphical solution of the eq. (1.2) is based on using the vectors gi ðnÞ ¼ si n þ s=i t, i¼c,f. Here s=i ðtÞ- is the angular derivative of free surface energy and t-is the unit vector, situated on the plane Pi touching the Wulff figure at a point A (Fig. 2), n-is the unit vector and is the normal to the plane Pi. It is convenient to rotate all vectors in eq. (1.2) by 90 , so this balance condition at the TPL described with the help of the vector equation [4,5,11]:   gS ðnS Þ ¼ gf nf þ gm ðnm Þ

ð1:3Þ

In Fig. 3 are shown an example of the graphical solution of this eqs. (1.3) according to CahnHoffman approach (a) and correspondingly the construction, proposed by Voronkov (c) for any equilibrium position of interface surfaces at the three-phase line (TPL) (b). Note, that it is permissible to represent the Wulff figures Gi, i¼c,f,m, only qualitatively since unknown are the values of the interphase energies and the shape of polar diagrams. However, it is possible to take into account some general postulates in the case, when the singular orientations are absent during the crystal shaping. Since sc>sf, the Wulff figure Gf for the crystal front surface always will be enclosed within the figure Gc. At the same time the figure Gm for the liquid phase is represented by a simple circle. The centre of this figure Gm situated on Gc or on Gf. The position of circle Gm centre is defined by choosing an azimuthal orientation of one of the crystalline surfaces. If the crystallization front orientation is fixed, the centre of Gm is positioned on Gf and according to the Voronkov rule the unique solution eq. (1.3) is determined by an intersection point of the figures Gm and Gc. Only the “entering” point into Gc (if the movement around the circle Gm is counter-clockwise) represents the correct solution [4,5]. Thus, as a result of “moving” the Gm centre around figure Gf, is it possible to construct the angular diagrams 4c(4m) and 4f(4m), which allow one to define the orientations of two surfaces by using the orientation of the third surface [4e7] (see Section 8.1). If the sectors on figures Gf and Gc, are situated far from singular orientations, then it is possible to approximate its shape by simple circles. It important, that the distance between Gf and Gc in this case is practically constant and it is evident, that the angle 3¼4m4c is almost constant. The value of this angle has been determined in articles [2,14] as the function of surface tensions. The solution of the equation set (1.2) for the surface forces and moments of these forces balances gives the expression:

Fig. 3. Graphical solution of the Herring eq. (1.2) (a)-according to Cahn-Hoffman [11] and (c) according to Voronkov [4,5]. (b)-is the equilibrium arrangement of the interfaces in the vicinity of the three-phase line L. In Voronkov’s graphical construction (c) the Wulff figures Gf and Gc have a common centre O. The vector gm is equal to the radius of the circle Gm, which is the section of the spherical Wulff figure for the isotropic melt surface. The Cahn-Hoffman approach differs from the Voronkov approach also in the choice of the normal direction to the growth front and to the melt interface (see [5]).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

 3 ¼ arccos

s2m þ s2c  s2f 2sc sm

5



ð1:4Þ

=

=

Plus any correction term by taking into account the angular derivatives sc and sf . Quantitative estimations allow one to neglect these terms due to their small values (only if the crystalline surfaces at the TPL are ordinary and not singular facets [14]). (4) According to Bolling and Tiller [12] the growth front is always curved due to the Gibbs -Thomson effect in the vicinity of line L and does not coincide with an isothermal surface T0 (Fig. 4). According to Voronkov [3], the lateral crystalline surface can also be curved in the vicinity of the TPL like the growth front due to the surface mass transfer. Such Gibbs-Thomson distortion of the lateral surface in the vicinity of TPL is especially notable if the singular facet adjoins or is connected to the line L. This supercooling causes the generation and movement of elementary steps on such a singular facet (Fig. 4b). The singular facet lags behind the isotherm T0 where in the cold region the generation of two dimensional nucleation occurs. The difference of chemical potential between solid and the melt drives the atom stream from the melt along the lateral crystal surface. As a result the lateral crystalline surface is curved. The angle c0 between an initial position of lateral crystal surface at the TPL and its final position (far from L, Fig. 4b) is: c0 ¼ KV 1=3 dT

ð1:5Þ

Here dT- is the supercooling of the melt at the TPL, V- is the velocity of L, K- is a coefficient depending on the speed of surface self-diffusion. The “visual” angle between the tangent to the melt meniscus and the crystal surface far from the line L was called by Voronkov an “effective growth angle” 3* [3e7]. The value of this angle was determined as:

Fig. 4. The typical situation near the three-phase eline. The angles 4iare determined as deviations from the singular facet at the crystallization front. The growth front is concave and is intersected by a crystalline facet (on the right side). (a)-the vicinity of the three-phase-line L, far from the facet . Here the angle c0/0 . (b)-the vicinity of the facet with orientation n0. The growth front curvature at the line L must be in equilibrium in accordance with the Bolling-Tiller rule [12]. The crystalline steps move and accumulate near to the TPL. Note, in this case (b) the lateral crystal surface is curved. The angle c0 between the lateral surface at the TPL and the final surface position far from L is determined by eq. (1.5).

6

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

3 ¼ 4m  4C ¼ 4m  4C  c0 ¼ 3  c0

ð1:6Þ

In general, the angle 3* can be the negative, but the growth angle 3 between the tangent to the meniscus and the curved crystal lateral surface at the TPL (Fig. 4b) is always positive and varies very weakly due to the small magnitude of the angular derivatives s0i of crystalline surfaces i¼c,f. If a singular facet adjoins the line L, so the angle between the tangent to the meniscus and the new crystal lateral surface (far from the TPL) is variable and is determined by eq. (1.6). The Voronkov approach, described above, provides an explanation of the radial instability of Si whiskers growing by the Vapor-Liquid-Solid method [15,16]. After the formation of a Si-Au liquid drop on a Si substrate growth begins in a vertical direction forming silicon whiskers. The base of the whiskers always has the shape of a solid “cone”. The start of radial oscillations correlate with their dependence on the whisker radial dimension and appears when the whisker diameter 2r is less than 1m. [15] (see Fig. 5a). In this case the angle c0 between the tangent to the lateral crystal surface at the TPL and the crystallization direction V (Fig. 5b) depends not only on the meridional, but also on the second, radial curvature R1(t): c0 ¼ KV 1=3 dT þ K1 V 1=3 R1 sm =sc  0:5K1 V 1=3 R1 Correspondingly the effective growth angle 3* (see eq. (1.6)) in this case depends simultaneously upon the supercooling at the crystallization front dT and also on the whisker radius r. Since the supercooling dT depends on the supersaturation in the Si-Au solution-melt drop and this supersuturation also depends on the melt drop radius R(r), so such a crystallization system becomes unstable at small whisker diameters. In [16] the VLS processes can be described as a second order dynamical system so that the diameter oscillations of the submicron Si whiskers can be explained by the appearance of limit cycles on the phase-plot. The formation of the base cone always finishes with the first whisker diameter oscillation. It is evident from Fig. 5a, that the oscillation does not die when the whisker has a submicron diameter 2r<1m. So the theory of crystal surface formation, developed in [1e7], provides a detailed description of the external crystal morphology during crystal pulling from the melt. In general, the lateral crystal surface orientation (Fig. 4b) is described by the equation:

Fig. 5. (a)-photo of Si submicron whiskers [15]. (b)- The oscillation of the Si whisker diameter is explained by positive feedback and by the appearance of a closed cycle on the phase-plane portrait of the second order. Such a dynamic system model of Vapor-Liquid-Solid crystallization method can be see in [16].

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

4C ¼ 4m  3 þ c0

7

ð1:7Þ

and depends not only on the values of the boundary free energies and the initial seed crystal orientation, but is also dependent on many other factors (the temperature gradients, coefficients of surface mass transfer, kinetic coefficients at the growth front, capillarity etc). Generally, the angle 3* (between the tangent to the meniscus and the lateral surface of the crystal far from the TPL) is variable and is determined by the eq. (1.6). But the angle 3 is determined by the free energies of the boundary surfaces. This angle varies very weakly in the case of nonsingular boundary surfaces near the TPL and can be considered as a fundamental constant for a crystalmelt system. The conditions near the TPL in the Cz process can be considered in the isotropic approach, since for the typical temperature gradients at the crystallization front the singular facets are localized on a small section of the crystal perimeter. The conditions for a constant growth angle 3 is especially important for mathematical modelling of the dynamics of the Cz method and for the design of crystal diameter feedback control. This section is based on work in reference [17]. 2. Low order model (LOM) of the Czochralski method This section develops a linearized equation set, which adequately describes the behavior and provides a quantitative correlation between the input and the output of this dynamical system considered as a mathematical model in automatic control theory [23]. Three fundamental physical laws e (i) the mass conservation law, (ii) the heat balance and (iii) the thermodynamic equilibrium condition at the three-phase line and as a consequence - the constant growth angle condition - are the basis of the mathematical modelling. The LOM is based on linearization of these fundamental physical laws. It is important to keep the major factors controlling the dynamics of the crystal diameter variations near its crystallization front during any proposed simplifications. The goal is to develop a reduced set of linear differential equations. This low order mathematical model predicts the behavior of the controlled object under various circumstances without carrying out expensive tests of the real physical system, and it simplifies the work needed to design a system for controlling the crystal diameter in Czochralski (Cz) and Liquid Encapsulated Czochralski (LEC). 2.1. Thermodynamic equilibrium condition at the three-phase line and capillarity The slope angle of the tangent to the crystal lateral surface near its crystallization front at the three-phase line b(t)¼a(t)30, is measured from a vertical line (Fig. 6), were a(t) is the slope angle of the tangent to the meniscus and 30 is the growth angle, which is the physical constant for the crystal material. The constant growth angle condition 30¼const allows one to write the well known equation: r_ ðtÞ ¼ Vc tgbðtÞ

ð1Þ

_  HðtÞ, _ Since the crystallization rate is equal to Vc ðtÞ ¼ V0  hðtÞ were V0 is the crystal _ and HðtÞ _ pulling rate, hðtÞ are the rate changes of meniscus height and melt level in a crucible, so linearization of eq. (1) with second order accuracy gives: d_r ðtÞ ¼ Arr ðtÞdrðtÞ þ Arh ðtÞdhðtÞ þ Arh_ ðtÞdh_ þ ArH_ ðtÞdH_ þ ArV ðtÞdV0 ðtÞ

ð2Þ

8

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 6. Scheme of the basic Cz method with notations used in the text. The meniscus and some parts of the growing crystal are allocated in the LOM as an observed region.

Here the coefficients Aij(t) are equal to: Arr ðtÞ ¼ Vc ðtÞb0r =cos2 b Arh ðtÞ ¼ Vc ðtÞb0h =cos2 b Arh_ ðtÞ ¼ ArH_ ðtÞ ¼ tgb ArV ðtÞ ¼ tgb

ð3Þ

The partial derivatives b0r and b0h can be determined by use of the meniscus height h(r,a), which depends on the current crystal radius and the tangent slope angle to the crystal lateral surface at the three-phase line, since b0h ¼ a0h ¼ 1=h0a and b0r ¼ a0r ¼ ðh0r =h0a Þ. Note that these expressions are correct for every condition of crystal pulling. There are several approximations of the capillary Laplace equation with reference to the Cz crystallization method (see Section 8 of this review). The meniscus shape and height are determined by the diameter 2r(t) of the three-phase line and the melt surface tension sLG only and are given by the well known Tsivinskii approximation [18]: 1=2  ð4Þ hðr; aÞ ¼ a2 ð1  sinaÞ þ a4 cos2 a=16r02 a2 cosa=4r0 After differentiation the required partial derivatives are:  a2 sina a cosa 2r cosa 1 B B h0a ¼ 4r 4r a sina i h a2 cosa a cosa h0r ¼ B 1 2 4r 4r

ð5Þ ð6Þ

where B¼(1sinaþa2 cos2a/16r2)1/2. This approximation more precisely represents the meniscus height when crystal radius r(t)>a exceeds the melt capillary constant a¼(2sLG/rlg)1/2 [19]. In Fig. 7 are shown the menisci, calculated by use of the Tsivinskii approximation.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

9

Fig. 7. Tsivinskii approximation. (a) meniscus heights, calculated by eq. (4) for the growth angle 30¼10 , (b) the menisci profiles and (c) the tangent slope angle to the meniscus q(hi/a) with e0¼10 and r0/a¼4, (d) the tangent slope angle to the meniscus q(ri/a).The calculations used eqs. (7) and (9). Note that in eqs. (7) and (8) Hurle defines the capillary constant as a¼(sLG/rlg)1/2.

Eqs. (4)e(6) are generally correct, i.e. for pulling by the Cz method involving an initial, and final cone together with the cylindrical crystal rod. Under stationary conditions of crystallization involving a cylindrical crystal r¼r0¼const, h¼h0¼const, r_ ¼ h_ ¼ 0 and b(t)¼0, so Vc ¼ V0  H_ o and the eq. (3) simplified. For the control methods which use light reflection from the menisci it is desirable to know the relationship defining the menisci shape in the Cz method. The equation approximates well

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

10

to the meniscus profile curves obtained by Hurle [20] which were based on the Tsivinskii approximation : x ¼Rþ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2=A  H 2  2=A  y2  pffiffiffiffiffiffi lnfFðyÞg 2A

ð7Þ

Here we note x¼r/a, y¼h/a,R¼r0/a, H¼h0/a, a0¼90 30, and A ¼ 0:5ð1 þ BÞ; B ¼ sinða0 Þ=RH; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi H ¼ 2ð1  cosða0 ÞÞ þ B2  B;   pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 !# pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y 2 þ 2=A  H 1=2 2 2 pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln abs FðyÞ ¼ R þ 2=A  H  2=A  y  ð2AÞ H 2 þ 2=A  y2

ð8Þ

In Fig. 7b are shown the menisci, calculated using eq. (7) for a growth angle 30¼10 in nondimensional coordinates x¼r/a and x¼y/a. By use of the equation A¼0.5(1þB)in (8) it is possible to calculate the tangential slope to the meniscus against to its height q(hi/a) as:   ð9Þ q ¼ arccos 1  Ay2 Or calculate the angle q(ri/a) as function of the meniscus cross section radius at the corresponding height hi/a (Fig. 7).

2.2. Analysis of mass balance in the basic Czochralski method Let’s denote M(t), m(t) and m(t) the current values of the melt mass in a crucible, the mass of growing crystal and the corresponding meniscus mass. At any time the total mass in the basic Cz method is constant M(t)þm(t)þm(t)¼const. Thus from the total mass variations dM(t)þdm (t)þdm(t)¼const follows the relationship: _ _ þ dmðtÞ dmðtÞ _ ¼ dMðtÞ

ð10Þ

The mass of a crystal pulled from the melt upto the moment t can be represented by the expression: Z t mðtÞ ¼ prs r 2 ðtÞVc ðtÞdt  rs UðrðtÞÞ ð11Þ 0

Here U(r(t)) represents the volume of the crystallization front under the three-phase line. In the case of a flat crystallization front U(r(t))¼0, when it is convex this volume is positive þU (r(t))s0 in eq. (11) and is negative U(r(t))s0 for a concave crystallization front. In general, the meniscus mass is determined by the equation: mðtÞ ¼ rl ½Vm ðtÞHUðrðtÞÞ ¼ prl r 2 ðtÞhðtÞ þ prl a2 rðtÞcosaðtÞHrl UðrðtÞÞ

ð12Þ

Here Vm - is the meniscus volume in case of a flat crystallization front, a is the capillary constant and a(t)¼b(t)þ30 (see Fig. 6). The relationship Vm¼prlr2hþprlra2cos(a) was proposed earlier, see, for example, the fine survey of capillarity problems in the work of

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

11

Borshforth [21]. Consider for simplification the simple cylinder crucible shape with an inner radius R¼const. After linearization of the eqs. (11) and (12) it is possible to obtain: _ _ dMðtÞ ¼ prl R2 dHðtÞ

ð13Þ

_ dmðtÞ ¼ prs r½2Vc drðtÞ þ rdVc ðtÞ  rs U0r d_r ðtÞ

ð14Þ

0 _ r ðtÞ dmðtÞ _ ¼ m0r d_r ðtÞ þ m0h dhðtÞHr l Ur d_

ð15Þ

The derivatives in eq. (15) are defined by use of the eqs. (5) and (6) as:   m0r ¼ prl 2rh þ a2 cosa  a2 ra0r sina   m0h ¼ prl r2  a2 ra0h sina

ð16Þ ð17Þ

After substituting the relationships (13)e(17) into eq. (10) we obtain the second equation of the linearized model of the basic Cz method as: _ þ AHV ðtÞdVc ðtÞ _ dHðtÞ ¼ AHr ðtÞdrðtÞ þ AHr_ ðtÞd_r ðtÞ þ AHh_ ðtÞdhðtÞ

ð18Þ

The coefficients in eq. (18) are equal to: AHr ðtÞ ¼ 2rs Vc r=D AH_r ðtÞ ¼ ðm0r Hðrs  rl ÞU0r Þ=pD AHh_ ðtÞ ¼ ðprs r 2  2m0h Þ=pD AHV ðtÞ ¼ prs r 2 =D

ð19Þ

Here D ¼ rl R2  rs r 2 . The problem of the shape of the crystallization front when modelling the Liquid Encapsulated Czochralski (LEC) method was considered in the work of Johansen [22]. His approach and conclusions can be completely transferred to the common Cz method [23]. Note, that the obtained relationship in eqs. (14) and (15) allows one to take into account any “composite” shape of a crystallization front, which consists of convex and concave regions simultaneously. It is possible to take into account any different assumptions concerning the behavior of the volume U(r(t)) at the time of crystallization, which can be based on a posteriori analysis of crystal sections. In practice the change in shape of the growth front, causes negligibly small variations to its volume so it is possible to neglect this effect and accept in calculations U’ (r(t))¼0, as for example, for Lithium Niobate (LiNbO3) crystals. More importantly this effect can be used in the case of the growth of semiconductor crystals (GaAs, GaP, InP) by the LEC method, especially at the time of the initial stage (seeding stage). Since the obtained relationships (11), (12) and (19) allow one to take into account any arbitrary composite shape of a crystallising growth surface it is sufficient to add any realistic model U’(r(t))s0 in the calculations of the program task and to transfer processes in the open and closed states of a dynamical system. 2.3. Analysis of the thermal conditions Heat transfer analysis at crystallization is reviewed in numerous publications; a detailed analysis is not presented here. A detailed knowledge of temperature fields is required, first of all, for the calculation of thermoelastic tensions and for assessing the point defects and

12

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

dislocations in a crystal. The distribution of impurities in a crystal and the shape of the growth surface are determined by the hydrodynamic melt flows in a crucible. These streams are fundamentally important for the calculation of impurity transport to the growth surface and its distribution in the crystal. At present it is very difficult, or even impossible, to control convective flows and diffusion processes in real time during Cz crystal growth. But it is well known, that a constant crystal diameter is related to the stationary growth conditions near the crystallization front. It is necessary to maintain those conditions which produce the best crystal quality, but such growth conditions are only obtained experimentally. It is necessary to use averaged heat flows through the boundaries when modelling the Cz, LEC and Floating Zone methods. The variations of these streams correlate with the crystal diameter and melt meniscus changes. It is fundamentally important, that the crystal diameter variation is localized near the crystallization front. As a result of the localization perturbations it is impossible to use a simple linearization over coordinate r of the undisturbed thermal task solution. In order to find the streams at the crystal and meniscus boundaries using traditional methods it is necessary to determine the values of temperature at each point in space and time. Such solutions may now be obtained numerically. Unfortunately, the numerical methods that give excellent results in the heat- and mass transfer tasks, do not fit the task of real time control and are not used yet in the engineering practice of crystal growth. But the “photographic” map of a temperature field calculated using special software, for example the Matlab program package, can help in understanding and in verifying the modelling approach. Thus, it is necessary to derive a low order mathematical model of the thermal problem, appropriate for the Cz crystallization method. The simplifications are based on averaging methods and the asymptotic expansion of mathematical physics equations over small parameters. For this the following approximations are used in the LOM: 1. Not taking into account the presence of impurities and convective flow, for example, the Marangoni flow. Heat transport in a crystal, having an axial symmetry occurs only via thermal conductivity. A crystal is pulled in the z direction at a constant rate V0, the thermal conductivity coefficients, the thermal capacity and the crystal and melt densities are constant as also are the heat-transfer coefficients of the isotropic surfaces. Also not discussed are the problems of fluid and gas dynamics caused by the free and forced convections due to crystal rotation. The crystal rotation rate is taken to be constant. 2. An isotropic approach is used to define interface boundaries. This means, that crystal surfaces faceting is not taken into account. The crystal growth surface follows the shape of the crystallization isotherm due to the normal crystallization mechanism, but with the exception of a local region near to the three-phase line (Fig. 6). The size of such a local deflected region does not exceed 10O100mk and its presence is explained by the GibbsThomson effect [3]. Also not taken into account is the local faceting of a lateral crystal surface as a result of a local singular crystalline surface forming at the growth surface adjoining a three-phase line. So the macroscopic crystallizing growth surface curvature follows the shape of the isotherm Tc¼const and adjusts to the slow changes of the thermal fields in the growth zone, including the environment of the growing crystal. The microscopic curvatures of the solid surfaces near to the three-phase line always maintain their equilibrium status due to the fast crystallization kinetics. 3. An important simplification in the heat problem statement is the assumption about the small curvature of the crystallizing growth surface, which allows temperature averaging in a crystal to occur.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

13

4. The thermal conductivity in a liquid phase near to the solidified front does not exceed the thermal conductivity in a solid phase. The values of thermal conductivity of different crystalline materials at the temperatures near to the melting temperature are measured and well-known, but for the melts it not so. The presence of the convective mixing of melts in experiments increases its measured thermal conductivity compared with the thermal conductivity of corresponding crystals. In the proposed model of crystallization processes the melt region near to a solidified front in which the contribution of the convective mixing to thermal conductivity of a liquid phase can be neglected. A good estimation of thermal conductivity coefficients of melts in this model is the simple extrapolation curve of known datas for solid crystals in the field of temperatures which is a little higher than the melting temperature. The violation of this rule in calculations brings incorrect results in the physical point of view. 2.3.1. The general problem: definition In the classical determination of the Stephan-Boltzmann task the heat balance condition at the crystallization front Z¼S(r,t) is: ls Vn Ts ¼ ll Vn Tl þ rs Vc L

ð20Þ

Here ls, rs, ll,rl are the thermal-conductivity coefficients and the densities of the crystal and melt, VnTi is the corresponding temperature gradient in the line of the external normal to the growth front, Vc the crystallization rate, and L is the latent heat of phase transition. This equation is correct at any time. Taking into account a relation for differential along the external normal to a boundary v=vn ¼ ðv=vrÞcosðnrÞHðv=vzÞcosðnzÞ it is possible to rewrite eq. (20) as: ls ðS0r jTsr0 þ jTsz0 Þ þ ll ðS0r jTlr0 þ jTlz0 Þ ¼ rs Vc L

ð21Þ

1=2 where j ¼ ð1 þ S02 , Tir0 ¼ ðvTi =vrÞ j Tiz0 ¼ ðvTi =vzÞ for i¼s,l. Since in the approach of r Þ a nearly flat crystallization front S0r /0, from eq. (20) the evident conclusion follows where the main contribution to the heat balance introduced by the axial temperature gradients in the solid Gs ¼ Tsz0 and liquid Gl ¼ Tlz0 phases. For small deviations from a stationary crystallization state the crystal radius changes by a small value dr and the corresponding meniscus height to dh. The new temperature gradients at the crystallization front are Gs and Gl and with second order accuracy can be represented as:

Gi ¼ Gi0 þ G0ir dr þ G0il dh þ 0ð3Þ i ¼ s; l

ð22Þ

Here Gi0 are the gradients in an initial stationary state. Here it is possible represent the current crystallization rate as: Vc ¼ Vc0 þ dVc

ð23Þ

Initially, for steady crystal pulling with a constant rate V0 without any changes of its cross section the crystallization rate is equal to Vc0 ¼ V0  H_ 0 . Since the meniscus height as well does not change, so h_ ¼ 0 and from the mass balance condition the well-known expression H_ 0 ¼ rs r02 V0 =ðrl R2  rs r02 Þ follows for the simplest case of a right cylindrical crucible with constant radius R¼const. But in the general case it is necessary to take into account the meniscus height variation. Since the crystallization rate is determined by the expression:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

14

Vc ¼ V0 

dSðr; tÞ ¼ V0  S0t  S0r rt0 dt

ð24Þ

and S(r,t)¼H(t)þh(r,t) so taking into account the above simplification S0r /0 we get the _ _  hðtÞ and consequently: crystallization rate in the general case as Vc ðtÞ ¼ V0  HðtÞ dVc ¼ dV0  dH_  dh_

ð25Þ

Note, the co-ordinates of the point of origin are connected with the bottom of the moving crucible (Fig. 6). After substituting (22) and (25) into the heat balance eq. (21) we find the relation which connects the crystal radius and meniscus height variations with the variations of the crystal rate: Ir dr þ Il dh ¼ rs LdVc

ð26Þ

Here Ij ¼ ll G0lj  ls G0sj , j¼r,h. For the analysis of transfer processes in the controlled object, that is the investigation of the dependences dr(t) and dh(t), one expects the determination of explicit functions G0i in (26). Following perturbations of steady crystallization conditions, for example after a change of heating power or crystal pulling rate, there appears a new crystal area having its length l0 that is different from the right cylinder (see Fig. 6). Using this length l0 as the scale unit in the region of interest, the rest of the crystal will be considered as removed infinitely far from the crystallization front. The heat transfer solution in the crystal and the melt in a local region near to the crystallization front (shown in Fig. 6 as the observed area) is required. Non-dimensional variables x¼r/r0, y¼z/l0, z¼z/h0, v¼Vc/V0 and the non-dimensional temperature Ui¼(TiTa)/(TcTa), i¼s,lare used for the solid and liquid phases. Here r0- is an initial cylindrical crystal radius and h0-the meniscus height for steady state crystal pulling with this radius and constant pulling rate V0, Tc is the constant crystallization temperature and Ta is the surrounding temperature for the meniscus and, correspondingly, the small crystal part, having length l0. So eq. (20) for heat balance at the crystallization rate can be rewritten in a non-dimensional form as: Usl0  ðll l0 =ls h0 ÞUlz0 ¼ tðV0 =l0 Þv

ð27Þ

In eq. (27) t ¼ l20 L=cs CðTc  Ta Þ is an important variable, which describes the time scale of a transfer process (growth of crystal region l0), were C is the heat capacity and c¼ls/rsC is the crystal diffusivity. Now it is possible to determine the explicit expression G0ij in eq. (22). With this purpose it is necessary to solve the thermal conductivity equation: ðv=vt þ Vc v=vzÞTi ¼ ci divðgradTi Þ

ð28Þ

for each phase i¼s,l. By the use of dimensionless time t*¼t/ti, i¼s,l and introducing above the other dimensionless variables, we rewrite this equation for the solid phase, having the coordinate’s origin z¼0 at the crystallization front, as:   l20 vus l20 V0 vus v2 u s 2 1 v vus ¼ 3 þ 2 x þ v s  cs ts vts cs l0 vy vy x vx vx

ð29Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

15

and for the liquid phase accordingly:   h20 vul h20 V0 vul v2 u s 2 1 v vul ¼ 3 þ 2 þ v x l  cl tl vtl cl h0 vz x vx vx vz

ð30Þ

Here tsand tl are the reference times, scaling duration of thermal processes in a solid and liquid phases and 3s¼l0/r0 and 3l¼h0/r0 -are the small parameters which have appeared in the task. In a sense in our problem it is necessary to observe the thermal processes in the crystal and melt regions in the time scale of the crystallization transfer process duration t. In other words in eqs. (29) and (30) consider, thatts¼tl¼t. This simple estimation shows, that in this time scale the thermal processes can be observed as a quasi-stationary approach. The estimation of the contribution of the time dependent component should be defined by comparison with the values of the components in right and left-hand parts of these eqs. (10) and (11). Since l20 =cs t ¼ CðTc  Ta Þ=L ¼ 3Ts and h20 =cl t ¼ 3Ts h0 cs =l0 cl ¼ 3Tl it is evident, that the stationary condition, for example in the solid phase, is fulfilled at 3Ts << 32s or otherwise l0 >> ð3Ts r0 Þ1=2 . An example using known data: for a Germanium crystal 3Ts¼ 0.00123*DT (Ge data: Cp¼0.74 cal g1K 1, rsL¼ 601 cal cm3 and rs¼ 5.26 g cm3 get), for Gallium Arsenide it is even less 3Ts¼ 0.000576*DT (the GaAs data: Cp¼0.4184 Jg1K1, L¼ 726.76J g1) and for Gadolinium Gallium 3Ts¼ 0.00129*DT Garnet (the GGG data: Cp¼ 0.0586 J g1K1, L¼ 455.4 J g1). Taking into account, that under real growing conditions the temperature of the surroundings in a local region (Fig. 6) does not much differ from the crystallization temperature at the front, we shall accept that Ta/ Tc0.9. Ta/Tc0.9. If the absolute value DT100  C, so for a crystal with r0z(10O50)mm the temperature distribution is formed and fixed in time for an additional crystal growth of length l0z(1O2)mm. The conditions of the quasi-stationary state in the melt under the crystallization front are also fulfilled 3Tl3Tsand h0>>(3Tlr0)1/2 since for the crystal radius r0z(10O50)mm the meniscus height h0z1.2a and at the typical values of capillary constant h0z(4O6)mm. The convective component of a thermal flow in eqs. (29) and (30) is determined by the Peclet number Pe ¼ r0 VC =ci, i¼s,l and depends on the crystallization rate VC(t), i.e. from the displacement of the meniscus melt to the crystallization front. In the first approximation consider, that the melt displacement rate to the crystallization front is equal to the crystallization rate Vc. Strictly speaking, in general, the value of the parameter l20 V0 v=cs l0 ¼ 3s Pes x3l Pel is variable and it depends on the moment in time of the transfer process. Since the crystallization rate can accept both values VC>V0 and also VC
G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

16

By use of the dimensionless variables the thermal task eq. (29) can be written as:   v2 Us vUs 2 1 v vUs ¼0  3s g þ 3s x vy2 vy x vx vx with the corresponding boundary conditions:

ð31Þ

us jy¼0 ¼ 1; 0  x  rð0Þ

ð32Þ

u0sn þ Bius jrðyÞ ¼ 0

ð33Þ

Here u0sn is the temperature gradient along the outer normal line to the lateral crystal surface, Bi ¼ ar0 =ls is the Biot number, a- is the effective coefficient of heat exchange from a lateral crystal surface and its surroundings. The eq. (32) specifies the constant temperature at the crystallization front and eq. (33) determines the heat exchange from the lateral crystal surface. The intensity of this heat exchange is determined by the Biot number. Note here, that in our analysis the boundary condition at the upper crystal end ( y/N) is not essential because the problem consists in the determination of a gradient derivative function G0sj , j¼r,z, instead of the crystal temperature itself. Assume that the front curvature is small. This allows one to use the averaging temperature us(x,y) in every cross-section S( y)¼pr2( y). Introducing a new R rðyÞ variable TðyÞ ¼ 2p 0 us ðx; yÞxdx=SðyÞ. Assuming us ðx; yÞjrðyÞ ¼ TðyÞ, after term-wise integration of eq. (31) becomes: Z Z

rðyÞ

ð34Þ

u00sy sy xdx ¼ 0:5r 2 ðyÞTy00 ðyÞ

ð35Þ

0 rðyÞ 0

Z

u0sy xdx ¼ 0:5r 2 ðyÞTy0 ðyÞ

rðyÞ  0

  rðyÞ 1 v vus vus vus xdx ¼ x ¼ rðyÞ jx¼rðyÞ x vx 0 vx x vx vx

ð36Þ

Taking into consideration that the slope of the lateral crystal surface is determined by an angle b(r) (Fig. 6), so the tangent on a three-phase line is equal to rz0 z ¼ tgb and the normal

Fig. 8. The local part of a growing crystal (separated from the melt). See details in the text.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

17

derivatives to it is v=vn ¼ cosbðv=vrÞ  sinbðv=vzÞ. Since it is possible to rewrite the boundary conditions (33) as: 0 0 u0sx ¼ Bius þ 32 s ry ðyÞusy ; x ¼ rðyÞ

ð37Þ

By use of the eqs. (34)e(36) one can obtain for the new “disturbed” crystal section instead of eq. (31) a new one-dimensional heat equation:   v 2 vTðyÞ vTðyÞ ð38Þ r ðyÞ  3s gr 2 ðyÞ  2Bi32s rðyÞTðyÞ ¼ 0 vr vy vy Here Bi ¼ ar0 =ls cosb. Now eq. (38) and the corresponding boundary conditions describe a crystal quasi-stationary state, which is formed in time t in the crystal section l0. eq. (38) is fulfilled for every crystal cross-section of this new disturbed crystal section. Our interest consists in the determination of the temperature gradient variation dGs as a function of any small variation dr of the crystal radius. After linearization of eq. (38) with regard to the parameter variations Pes ðtÞ ¼ g0 þ dgðtÞ due to the variable crystallization rate Vc ðtÞ ¼ Vc0 þ dVc ðtÞ, one can get the following equation:   2 0 0 2 2Bi 2 Bi dGsy  3s g0 dGs þ 3s dT ¼  ðGs0 drÞy þ 3s g0 Gs0 dr þ 3s T0 dr ð39Þ r0 r0 r0 with the boundary condition: dTð0Þ ¼ dTðNÞ ¼ 0

ð40Þ

This last eq. (40) determines the condition of the temperature persistence at the crystallization front and at a remote distance from it ( y>>l0). We can note here, that during the transfer processes what is more significant is the effect of the convective heat exchange with the crystal environment, i.e. the component 3s g0 dGs , as compared with the heat exchange through the crystal lateral surface, i.e. ð32s 2Bi=r0 ÞdT in eq. (39). This fact results from, first of all, the locality of crystal surface variation. Now let’s solve this equation (39). To a first approximation, we shall not take into account the components with small parameters 3s/0. In this case eq. (39) become simpler, namely it reduces to dG0sy ¼ 2ðGs0 dr=r0 Þ0y and its solution looks like: dGs ð0Þ ¼ ð2Gs0 =r0 Þdrð0Þ

ð41Þ

From this relationship it follows, that G0sr ð0Þ ¼ 2Gs0 =r0 . Now we shall keep in eq. (39) the components of first order of smallness 3s. But again the dependence, described by eq. (41) is the solution of the derived “extended” equation dG0sy  3s g0 dGs ¼ r20 ððGs0 drÞ0y þ 3s g0 Gs0 drÞ under the boundary conditions eq. (40). So with second order accuracy the variations of an average axial temperature gradient in the solid phase at the crystallization front are independent of the heat exchange conditions of a crystal with its environment and can be presented by the dependence (41). The value Gs0 is the initial temperature gradient at the crystallization front before the crystal diameter variation. The derived equation (41) completely coincides with a relationship obtained earlier in the work of Voronkov [24]. In his work assumptions were used about the thermal environment invariability around the crystal with its local variation in diameter and as a consequence - the

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

18

constancy of the heat flow through the crystallization front, i.e. pr02 ðtÞGs ðtÞ ¼ const. A simple linearization of the last expression reduces to the above relationship (41). Consideration of the thermal analysis in a crystal includes the heat exchange from a lateral surface and the exchange due to displacement of a solid phase. It was found, that the heat exchange through its lateral surface should be taken into consideration only at a large Biot number Bi ¼ ar0 =ls cosb in eq. (39). In other words, with very intensive heat exchange due to a local crystal shape variation. For the remaining cases (really for all practically significant growing conditions) the relationship (39) as represented is correct. Note, the basic propositions in the analysis were the assumptions about the small macroscopic curvature of the crystallization front. 2.3.3. Analysis of the thermal problem for the meniscus In order to determine the values and resulting form of the derivatives G0lr and G0lh in eq. (22) it is necessary to take into account the geometrical features of the menisci, rising over the free melt surface and up to the growing crystal. This melt meniscus can be considered as the thin round disk (see Fig. 9) which is bounded at the top by the nearly flat crystallization front and from below by the free melt surface (the meniscus base). To a first approximation, it is possible to neglect the curvature of a meniscus lateral surface since r(t)>h(t) at the initial growth stage (seeding) where the crystal radius r0>>h0w1.2a under steady state growth of a massive crystal (see Fig. p 7).ffiffiffi The melt capillary constant a ¼ ð2sLG =rl gÞ1=2 has typical values a¼2O4mm and h0 ¼ 2a is the asymptotic meniscus height in the Cz method at 30¼0 . It is also possible to carry out the analysis of the thermal task in a liquid phase using asymptotic decomposition (see the work Zino, Tropp [25]). Let’s keep the convective component of a thermal flow in eq. (30) in a similar manner to the above and represent the parameter 3lPel¼3lg, were g¼O(3):   v2 u l vul 2 1 v vul þ 3 ¼0 ð42Þ x  3 g l l vz x vx vx vz2 The boundary conditions in this case are: ul jz¼1 ¼ 1

ð43Þ

ul jz¼0 ¼ ðT0  Ta Þ=ðTc  Ta Þ ¼ u0

ð44Þ

u0lx þ Biul jx¼1 ¼ 0

ð45Þ

Fig. 9. Liquid meniscus joining the “disturbed” crystal. Details in the text.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

19

Here T0 is the temperature at the meniscus base, Ta is the environment temperature near to meniscus, Bi ¼ ar0 =ll is the Biot number, which depends on the effective heat exchange coefficient a of the meniscus and its surroundings. The boundary conditions (43) and (44) determine the non-dimensional temperature at the upper and lower meniscus boundaries. The value of the Biot number depends on the heat exchange intensity at the meniscus lateral surface. Initially taking this coefficient as not small Bi¼0(1), in other words the heat exchange intensity is considerable. By use of an asymptotic expansion we represent the solution of task (42)e(45) as an expansion in a series [25]: uðx; zÞ ¼ u0 ðx; zÞ þ 3l u1 ðx; zÞ þ .

ð46Þ

After substituting eq. (46) into an initial eq. (42) we get a so called “external” solution, involving one main item u0 ¼ u0 ðzÞ where one may find as the solution of the homogeneous equation with corresponding boundary conditions: v2 u0 ðzÞ ¼ 0; u0 ð0Þ ¼ u0 ; u0 ð1Þ ¼ 1 vz2

ð47Þ

The remaining components of decomposition (46) are equal to zero. So in the central (core) meniscus part, the temperature distribution is similar to a linear distribution according to the dependence: u0 ðzÞ ¼ ð1  u0 Þz þ u0

ð48Þ

Now let us determine the contribution of the meniscus temperature distribution of the heat exchange from its surface with the intensityBi¼0(1). It is evident, that near to the lateral meniscus surface there is situated some boundary layer. For the determination of the temperature distributions in this layer it is necessary to introduce new “extended” coordinates near to the meniscus edge using the x coordinate transformation r ¼ ð1  xÞ=3l . At the same time the scale along the vertical axis z remains invariable. The “internal” solution for this boundary layer again involves searching as an expansion in the series: Wðr; zÞ ¼ 3l W0 ðr; zÞ þ 32l W1 ðr; zÞ þ .

ð49Þ

Now the investigated area has size limits 0
ð50Þ

ð51Þ

By taking into consideration the relationship ð1  3l rÞ1 ¼ 1  3l r þ oð32 Þ we get for W0 (r,z) from the task:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

20

v2 W 0 v2 W 0 þ ¼0 vr2 vz2

ð52Þ

W0 jz¼1 ¼ W0 jz¼0 ¼ 0; W0 /0 PpJ r/N

ð53Þ

W 0 j r ¼ 0 ¼ Biu0 ðzÞ

ð54Þ

z¼1 Note, that the quasi one-dimensional temperature distribution in a meniscus is obtained only if the Biot criterion, which defines the heat exchange surface with an environment, has the order O (32). With an intensive heat exchange 3l Bi ¼ 3l ar0 ¼ ah0 ¼ Oð1Þ, then the third item in the eq. (50) can be represented as Biðu0 ðzÞ þ 3l W0 ðr; zÞ þ oð32 ÞÞ ¼ Biu0 ðzÞ þ ah0 W0 ðr; zÞ þ oð32 Þ and the boundary condition (54) becomes: 0  ah0 W0 ¼ ah0 u0 W0r

ð55Þ

More exactly the temperature distribution in the boundary layer can be obtained from a solution ofPthe tasks (52)e(55) The general solution form is pffiffiffiffiffi by use pffiffiffiffiffiof the Fourier method. 2 C expð l ¼ ðnpÞ ; n ¼ 0; 1; 2. the constants Cn rÞsinð l zÞ; l W0 ðr; zÞ ¼ n¼N n n n n n¼0 may be obtained from the expansion of the Fourier series and by use of boundary conditions (54) or (55). So the temperature distribution in a meniscus ul ðx; zÞ ¼ u0 ðzÞ þ 3l W0 ðr; zÞ þ 0ð32 Þ is found as a one-dimensional temperature distribution having a small correction, which is essential only near to the meniscus lateral surface under intensive heat exchange from this surface. The width of the boundary layer may be estimated as p o(3 ffiffiffi l) or in dimensional form as o(h0), in other words comparable with the meniscus height h  2a, were a the capillary constant is specified as its upper bound. From the above analysis it follows, that the convective component of heat exchange appears in asymptotic solution form (49) only in items of second order smallness. The Marangoni convection can increase the boundary width in general, essentially more at the initial (seeding) stage of crystal growth, but even in this case our conclusions do not change for “large” crystals (r0>>a). Thus, with the approach of a flat crystallization front and not too much intensity of meniscus heat exchange from surface and its environment, the main temperature distribution (48) is linear along the meniscus height. And if the melt temperature at the meniscus base is equal to any T (0)>Tc¼const, so the axial temperature gradient in it is determined as: Gl ðr; zÞ ¼ ðTc  Tð0ÞÞ=h0 þ oð3l Þ

ð56Þ

Correspondingly, the temperature gradient derivatives in eq. (22) are: G0lr ¼ 0 G0lr ¼ ðTc  Tð0ÞÞ=h20

ð57Þ

In the Liquid Encapsulated Czochralski (LEC) method the melt surface is covered by a liquid encapsulant and consequently the Marangoni convection is vastly suppressed. The known experimental distribution of linear temperature in the center part of a melt meniscus is reported in many publications, in particular in the work Kuroda, Kozuke [26]. It is possible to get a more direct presentation about the temperature distribution in a crystal and an adjoining meniscus by use of the Matlab programmes for heat task solutions using the

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

21

Finite Element Method. In Fig. 10 is shown such a distribution, which was obtained using the software package pdetool of Matlab7. Finally, taking into consideration variations of the melt temperature at the meniscus base dT (t) and possible pulling rate variations dV0(t), after corresponding substitutions into Ij ¼ ll G0lj  ls G0sj , j¼r,h in the heat-balance eq. (26) we obtained the above derivatives (41) and (57) and eq. (25) for dVC(t) getting the third equation of the linear mathematical model for the Cz method as: _ ¼ Ahr ðtÞdrðtÞ þ Ahh ðtÞdh  dHðtÞ _ þ AhT ðtÞdTðtÞ þ dV0 ðtÞ dhðtÞ

ð58Þ

where the coefficients (called thermal coefficients) of the Low Order Model of the Cz method are: Ahr ðtÞ ¼ ð2lS GS0 ðtÞÞ=rS LrðtÞ Ahh ðtÞ ¼ ll ðTC  TðtÞÞ=rS Lh2 ðtÞ AhT ðtÞ ¼ ll =rS LhðtÞ

ð59Þ

It is necessary to note, that it is possible to extend the derived result to the case of nonstationary growth of rotationally symmetric crystals, in particular to the initial stage of cone growth. In general, the values of the solid phase gradient GS0 ðtÞ and the melt temperature are both slowly variable values during the growing process. Therefore eqs. (58) and (59) are written, in general, as time-dependent. At the same time the slowness of the change of heat conditions in a Cz growing zone allows the use of the “frozen” coefficients approach [27].

Fig. 10. The typical distribution of temperature is one-dimensional for a nearly flat crystallization front. Calculations were carried out using the finite element method obtained with the software package pdetool (Matlab7). The melt meniscus and the crystal after some transfer process are constructed as the two separated areas. Accordingly, the solidified front is represented by two horizontal boundary lines on which are noted the constant (crystallization) temperature. We used the Neuman boundary condition, the heat-transfer coefficient for lateral surfaces of a crystal and melt equal to 101.

22

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

So in the proposed Low Order Model (its thermal part) it is required to know or estimate only two variable values - the temperature gradient GS0 ðtÞ in the crystal at the crystallization front and the melt temperature overheating T(t)>Tc¼const at the meniscus base. 2.4. The standard form of the LOM for the Czochralski method (continuous time) The conservation laws (the mass and heat balance and the growth angle constancy) exist and are obtained for any crystallization time. The linearization of these three conservation laws yield three differential eqs. (2), (18) and (58), which form the equations set of the Low Order Model (LOM) for the Czochralski method: d_r ðtÞ ¼ Arr ðtÞdrðtÞ þ Arh ðtÞdhðtÞ þ Arh_ ðtÞdh_ þ ArH_ ðtÞdH_ þ ArV ðtÞdV0 ðtÞ _ ¼ Ahr ðtÞdrðtÞ þ Ahh ðtÞdh  dHðtÞ _ þ AhT ðtÞdTðtÞ þ dV0 ðtÞ dhðtÞ _ þ AHV ðtÞdV0 ðtÞ _ dHðtÞ ¼ AHr ðtÞdrðtÞ þ AHr_ ðtÞd_r ðtÞ þ AHh_ ðtÞdhðtÞ

ð60Þ

The coefficients Aij are presented in detail below. In the first equation for the crystal radius variation dr_(t) coefficients Aij are described as: Arr ðtÞ ¼ Vc ðtÞb0r =cos2 b ¼ A11 Arh ðtÞ ¼ Vc ðtÞb0h =cos2 b ¼ A12 Arh_ ðtÞ ¼ ArH_ ðtÞ ¼ tgb ¼ A13 ¼ A14 ArV ðtÞ ¼ tgb ¼ A15

ð61Þ

_ Respectively for the meniscus height variations dhðtÞ Ahr ðtÞ ¼ ð2lS GS0 ðtÞÞ=rS LrðtÞ ¼ A21 Ahh ðtÞ ¼ ll ðTC  TðtÞÞ=rS Lh2 ðtÞ ¼ A22 AhT ðtÞ ¼ ll =rS LhðtÞ ¼ A28

ð62Þ

and for the variation of melt level in a crucible dH_(t): AHr ðtÞ ¼ 2rs Vc r=D ¼ A31 AH_r ðtÞ ¼ ðm0r Hðrs  rl ÞU0r Þ=pD ¼ A34 AHh ðtÞ ¼ ðprs r 2  2m0h Þ=pD ¼ A35 AHV ðtÞ ¼ prs r 2 =D ¼ A37

ð63Þ

Here in eqs. (61)e(63) are shown also the reduced numerical notation of the coefficients Aij i¼1,.7, which are used in the numerical calculations . In general, the LOM is linear and time-dependent, first of all, due to the melt level changes in a crucible. For the basic independent variables of the LOM Cz method it is necessary to choose the vector XðtÞ ¼ ½drðtÞ; dhðtÞT ¼ ½x1 ðtÞ; x2 ðtÞT , since the third variable dH(t) derived from dr(t)anddh(t). It is evident, that for crystal pulling (by basic Cz) the melt level variation in a crucible is a consequence of crystal diameter and meniscus mass variations, not the contrary. Correspondingly the input vector in the LOM, in other words a vector of perturbations is UðtÞ ¼ ½dTðtÞ; dV0 ðtÞT ¼ ½u1 ðtÞ; u2 ðtÞT . In general, for the time dependent case, the matrix notation follows as:         a11 a12 b11 b12 x1 ðtÞ u1 ðtÞ x_ 1 ðtÞ ¼ þ ð64Þ a21 a22 b21 b22 x_ 2 ðtÞ x2 ðtÞ u2 ðtÞ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

23

For the stationary cylindrical crystal growth condition b(t)¼0 the model becomes autonomous (independent of time, certainly for limits of an assumption of the “frozen” coefficients) in the transfer process time scale:         a11 a12 0 0 x1 ðtÞ u1 ðtÞ x_ 1 ðtÞ ¼ þ ð65Þ a21 a22 b21 b22 x_ 2 ðtÞ x2 ðtÞ u2 ðtÞ In reduced form the LOM may be represented in the standard matrix form as: _ XðtÞ ¼ AXðtÞ þ BUðtÞ

ð66Þ

The matrix A(t)and B(t) and its components aijand bijare written in detail below in the general case of a time dependent model of the Cz method:     Arr þ Arh_ Ahr Arh þ Arh_ Ahh a11 a12 ¼ Arr AHr Arr AHr_ Ahr ArH_ AHr_ Ahh Arh AHr_ Ahh ArH_ AHr_ ð67Þ A¼ a21 a22 1þAHh_ 1þAHh_ 0 1 ArV ð1AHV Þ   A A rh hT 1ArH_ AHr_ b11 b12 A B¼ ¼ @ AhT ð1A _ AHr_ Þ ð68Þ rH 1AHV b21 b22 1þAHh_

1þAHh_

It is possible to define this linear equation set (66) as the model of the controlled object without inertia (the inertialess model) because the input melt temperature variation u1(t)¼dT (t) at the meniscus base does not take into consideration the real dynamics of melt temperature response due to changes of heating power. This question will be discussed later in Section 4. These melt temperature variations dT(t) represent some input system perturbation or temperature noise in the T-control channel in an open state. The pulling rate perturbations or controlled input in the V-channel u2(t)¼dV0(t)do not require consideration of any inertial properties, since they are direct and act immediately on a melt meniscus (see also Section 3). Let’s consider in detail weight control in the Cz method and determine the output equation for the Cz model in the case of an ideal weight sensor. It is possible to measure a growing crystal or a crucible with the melt. In the approach of an ideal weighing sensor both methods are completely equivalent. In detail for the real weight balance see Section 7. From the mass balance in the basic Cz method it follows that the time t is fulfilled as: Z t Z t _ ð69Þ FðtÞ ¼ prs r 2 ðtÞVo ðtÞdt þ mðtÞ ¼ prl R2 HðtÞdt o

o

It is obvious, that the rate (derivative) of measured mass can be the information about the current state of the independent variable vectors X(t). After differentiation and linearization eq. (69) with second order accuracy we get:     _ ¼ F_ o ðtÞ þ dF_ o ðtÞ ¼ prs r2 þ 2ro dr ðVco þ dVo Þ þ dm_ ¼ prl R2 H_ o þ dH_ ð70Þ FðtÞ o _ Since FðtÞ ¼ prs ro2 Vco ¼ prl R2 H_ o is the program rate of mass change, from eq. (70) it follows: _ ¼ dmðtÞ _ _ þ dmðtÞ dFðtÞ _ ¼ dMðtÞ

ð71Þ

24

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Replacing eqs. (13)e(15) into eq. (71) yields the output equation for weighing control: Y ¼ prl R2 dH_ ¼ c11 dr þ c12 dh þ d11 dT þ d12 dV

ð72Þ

Or in the standard matrix form the output equation in the Cz method (72) looks like: YðtÞ ¼ CXðtÞ þ DUðtÞ

ð73Þ

below are represented the matrices C(t)¼[c1,c2] and D(t)¼[d1,d2] for the time dependent Cz model: C ¼ ðc11 c12 Þ ¼ PK½ðAHr _ þ AH_ r_ a11 þ AH_ h_ a21 ÞjðAH_ r_ a12 þ AH_ h_ a22 Þ D ¼ ðd11 d12 Þ ¼ PK½ðAH_ r_ b11 þ AH_ h_ b21 ÞjðAH_ r_ b12 þ AH_ h_ b22 þ AHV _ Þ

ð74Þ

with PK¼prlR2. 2.5. Determination of the transfer functions for Cz (the LOM in continuous time) The Laplace transformations of eqs. (66) and (73) at zero initial conditions give: XðsÞ ¼ ðsI  AÞ1 BUðsÞ   1 YðsÞ ¼ CðIs  AÞ B þ D UðsÞ The explicit transfer functions can be obtain from the solution of the equation set: 0 1 10 1 0 ðs  Arr Þ ðArh þ sArh_ Þ sArH_ drðsÞ ArV dVðsÞ @ Ahr ðs  Ahh Þ s A@ dhðsÞ A ¼ @ AhT dTðsÞ þ dVðsÞ A sAHh_ s ðAHr þ sAHr_ Þ AHV dVðsÞ dHðsÞ

ð75Þ ð76Þ

ð77Þ

In detail the solution of these equations is shown below: WrT ðsÞ ¼ ½sArh_ AhT ð1 þ AHh_ Þ þ Arh AhT =DðsÞ WrV ðsÞ ¼ ½Arh ð1  AHV Þ þ Arh_ ðAhh  AHV Þ=DðsÞ WhT ðsÞ ¼ ½sðAhT þ Arh_ AHr_ Þ þ Arh_ AHr  Arr AhT =DðsÞ WhV ðsÞ ¼ ½sð1 n  AHV Þ  ðArr þ Arh_ Ahr Þð1  AHV Þ=DðsÞ

WHT ðsÞ ¼ s2 ½AHh_ AhT þ AhT Arh_ AHr_  þ sðArh AhT AHr_  Arr AhT AHh_ þ Arh_ AHr AhT Þ o. þArh AhT AHr sDðsÞ  WHV ðsÞ ¼ s2 ðAHh_ þ AHV þ 2Arh_ AHr_ Þ þ sðArh AHr_  Arr AHh_  Arr AHV  Ahh AHV Þ þ ArV ðAhr AHh_  Ahh Ah_r þ Ahr AHV Þ þ Arr Ahh AHV  Arh Ahr AHV  Arh_ Ahh AHr sDðsÞ þ Arh AHr

ð78Þ

Here D(s) - is the general LOM characteristic equation in the case of variable and time dependent coefficients: DðsÞ ¼ s2 ð1 þ AHh_ þ 2Arh_ AHr_ Þ  s½Arr þ Ahh þ Arr AHh_  Arh AHr_  2Arh_ AHr  Arh_ Ahr þ Arh_ ðAhr AHh_  Ahh AHr_ Þ þ Arr Ahh þ Arh AHr  Arr Ar_ H þ Arh Ahr  ArH_ Ahh AHr

ð79Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

25

In case of steady state cylindrical crystal growth r0¼const, b¼0, a¼30 j Arh_ ¼ ArH_ ¼ ArV ¼ 0. So the autonomous Cz model transfer functions can be determined from such an equations set: 1 10 0 1 0 Arh 0 drðsÞ ðs  Arr Þ 0 @ Ahr ðs  Ahh Þ s A@ dhðsÞ A ¼ @ AhT dTðsÞ þ dVðsÞ A ð80Þ ðAHr þ sAHr_ Þ sAHh_ s AHV dVðsÞ dHðsÞ Its solutions look like: WrT ðsÞ ¼ Arh AhT =DðsÞ WrV ðsÞ ¼ Arh ð1  AHV Þ=DðsÞ WhT ðsÞ ¼ sAhT  Arr AhT =DðsÞ WhV ðsÞ ¼ ½sð1  AHV Þ  Arr ð1  AHV Þ=DðsÞ 2 WHT ðsÞ ¼ fs n AHh_ AhT þ sðArh AhT AHr_  Arr AhT AHh_ Þ þ Arh AhT AHr g=sDðsÞ WHV ðsÞ ¼ s2 ðAHh_ þ AHV Þ þ sðArh AHr_  Arr AHh_  Arr AHV  Ahh AHV Þ o. þArr Ahh AHV  Arh Ahr AHV þ Arh AHr sDðsÞ

ð81Þ

The characteristic equation is simplified in this case and can be described as: DðsÞ ¼ s2 ð1 þ AHh_ Þ  sðArr þ Ahh þ Arr AHh_  Arh AHr_ Þ þ Arr Ahh þ Arh AHr þ Arh Ahr

ð82Þ

or in the general form as: DðsÞ ¼ s2  ða11 þ a22 Þs þ a11 a22  a12 a21 ¼ s2  ðspurAÞs þ detA

ð83Þ

2.6. LOM Czochralski in discrete time For the analysis and synthesis of a control system for the Czochralski (and LEC) method it is necessary define the discrete form of LOM. This transformation is based on the solution of the continuous model (66) in the Cauchy form Z t XðtÞ ¼ FðtÞXð0Þ þ Fðt  tÞBuðtÞdt ð84Þ 0

P where FðtÞ ¼ expðAtÞ ¼ ðAtÞv =n!; n ¼ 0; 1; .; N is the transfer matrix. n The analog-digital transformation of an input signal will be represented by a zero order predictor so the input and output discrete (with time T0) signals are extrapolated by step functions u(t). Since t¼(kþ1)T0, X(0)¼X(kT0)u(t)¼u(kT0) from eq. (84) we obtain: Z ðkþ1ÞT0 Xððk þ 1ÞT0 Þ ¼ FðT0 ÞXðkT0 Þ þ F½ðk þ 1ÞT0  tBdt ð85Þ kT0

Assuming the matrices A; B; C; D are invariable during the digitization time then after the introduction of new variables q ¼ ðk þ 1ÞT0  t, i.e. dq ¼ dt, q ¼ 0 at t ¼ ðk þ 1ÞT0 and q ¼ T0 at t ¼ kT0 it follows:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

26

Xðk þ 1Þ ¼ FðT0 ÞXðkÞ þ uðkÞ

hZ

0

i FðqÞdq B

ð86Þ

q

Or finally the discrete LOM form looks like the standard form:  b b Xðk þ 1Þ ¼ AXðkÞ þ BUðkÞ b b YðkÞ ¼ CXðkÞ þ DUðkÞ

ð87Þ

R b ¼ FðT0 Þ ¼ expðAT0 Þ; B b ¼ C and D b ¼ ½ T0 FðqÞdqB; C b ¼ D are the new where A 0 b it is matrices of this discrete model. For the determination of an explicit form of the matrix A b necessary to use the Hamilton e Kelly theorem [28], which gives in our case A ¼ a0 I þ a1 A. The constant multipliers ai can be obtained from the equations set:  expðl1 T0 Þ ¼ a0 þ a1 l1 ð88Þ expðl2 T0 Þ ¼ a0 þ a1 l2 Here l1,2- are the roots of the characteristic equation ðlI  AÞ ¼ 0. In the case of real root values the multipliers ai are equal to: a1 ¼ ½expðl2 T0 Þ  expðl1 T0 Þ=ðl2  l1 Þ a0 ¼ ½l2 expðl1 T0 Þ  l1 expðl2 T0 Þ=ðl2  l1 Þ b in the discrete model (87) can be determined as: So the matrix A     f11 f12 a1 a12 a0 þ a1 a11 b ¼ A¼ a0 a21 a0 þ a1 a22 f21 f22

ð89Þ

ð90Þ

where: l2  aii l1  aii fii ¼ expðl1 T0 Þ  expðl2 T0 Þ l2  l1 l2  l1 aij fij ¼ ½expðl2 T0 Þ  expðl1 T0 Þ here i; j ¼ 1; 2 and isj l2  l1

ð91Þ

In case of the complex root values l1;2 ¼ a  ib after proper transformations: fii ¼ expðaT0 Þ½sinðbT0 Þðaii  aÞ=b þ cosðbT0 Þ ð92Þ fij ¼ expðaT0 Þaii sinðbT0 Þ=b RT Since 0 0 ðexpðAtÞÞdt ¼ ½FðT0 Þ  IA1 , so the second matrix is described by the equation b  IÞA1 B. The matrix A1 B is equal to: b ¼ ðA B   1 a22 b11  a12 b21 1 EX1 EX2 a22 b12  a12 b22 1 A B¼ ¼ ð93Þ det a21 b11 þ a11 b21 a21 b12 þ a11 b22 det EX3 EX4 b components can be determined as: where det ¼ a11 a22  a12 a21 . So finally the matrix B     b ¼ ðdetAÞ EX1ðf11  1Þ þ EX3f12 EX2ðf11  1Þ þ EX4f12 ¼ B11 B12 ð94Þ B EX1ðf22  1Þ þ EX3f21 EX4ðf22  1Þ þ EX2f21 B21 B22

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

27

The quantization synchronism, i.e. neglecting a small delay in the work of a controller or computer, allows one to use for the discrete LOM the matrices C and D from the continuous time model. For simplicity below we do not use symbol L(cover) for the discrete model matrices. The basic mathematical model for the Czochralski method in its basic standard second order form allows one use modern automatic control theory for the analysis of the crystallization dynamics for the Cz method in an open and closed state. The LOM in this second order form uses a restricted number of parameters. Most are well-known and the pair parameters (the melt temperature overheating and the temperature gradient in solid phase at the crystallization front) can be estimated or defined experimentally. As an illustration in Fig. 11 are shown the graphics of LOM matrix A and B coefficients for a LiNbO3 crystal, calculated using eqs. (60)e(63) and parameters from Table 1. Also shown is

Fig. 11. Calculations of the LOM matrix coefficients for a LiNbO3 crystal using eqs. (60)e(63) and parameters from Table 1. (a) The calculated initial crystal cone and the selected crystal radius for which were calculated the coefficients aij(r) and bij(r) of matrix A and B. (b) For the selected crystal radius Ri ¼ 28 mm the values aij and bij are determined as the intersection of the horizontal line and curves aij(r) and bij(r).

28

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Table 1 Parameters for the LiNbO3 crystallization process Technological parameters (well-known) Seed crystal radius Crystal radius Radius of crucible Angle of crystal cone Crystal pulling rate

(m) (m) (m) (grad) (m/sec)

Rs R0 Rcr b V0

4.5*0.001 40*0.001 100*0.001 50 3*(0.001/3600)

LiNbO3 material parameters Thermal conductivity solid phase Thermal conductivity liquid phase Latent heat Growth angle Capillary constant Density of liquid phase Density of solid phase

(W/mK) (W/mK) (J/kg) (grad) (m) (kg/m3) (kg/m3)

ls ll L e0 A rl rs

1.1 0.6 178420 0 4/1000 3.8*0.001/(0.01)^3 4.65*0.001/(0.01)^3

Technological parameters (rough estimate for LiNbO3 crystallization process) Melt overheating (K) DT Temperature gradient in crystal (K/m) Gs

4 1500

the calculated initial crystal cone and the horizontal line, which corresponds to this initial cone section Ri ¼ 0.7r0 (28 mm) with the corresponding values aij(r) and bij(r). 3. Open loop analysis for the Czochralski dynamics By using eqs. (66) and (73) LOM in continuous time (or after z-transformation of LOM eq. (87) in discrete time) it is possible to represent the Cz method using the standard flow block, shown in Fig. 12: As will be shown later it is not difficult to calculate the crystallization rate and by measuring the sensor weight of the growing crystal its shape at any given pulling rate as well as the diameter of the cylindrical boule and crucible. For this purpose it is enough to know only two conservation laws: mass balance in the conserved system and the constancy of the growth angle. However is it possible to change and maintain in real processes these calculated velocities which depend on the thermal conditions in the growing zone. The LOM allows one to solve a problem of maintaining a state of established equilibrium in this zone or to carry out a calculated program task when the thermal conditions in the growth zone allow one to do it. So the LOM equations describe the behavior of the deviations XðtÞ ¼ ½drðtÞ; dhðtÞT ¼ ½x1 ðtÞ; x2 ðtÞT from any already established

Fig. 12. Vector-matrix diagram of the Cz method as the controlled object based on the LOM eqs. (66) and (73) in continuous time.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

29

stationary system state. The external actions on the system e the input perturbations in our case are the melt temperature and pulling rate variations UðtÞ ¼ ½dTðtÞ; dV0 ðtÞT ¼ ½u1 ðtÞ; u2 ðtÞT . The input vector U(t) in the LOM eqs. (66), (73) or in (87) at the same time include the input control for the temperature (T-) and rate (V-) channels when one designs a closed state of a controlled object. Formally, for the analysis of a controlled object in an open state (i.e. without feedback control) is not important to know the reason for the perturbations. What is significant is the fact of their presence and their mathematical description. For example, velocity perturbations u2 ¼ dV0(t) can be due to any puller drive defects or as a consequence of a stepwise slip of the seed crystal. Or, due to a non cylindrical crucible shape R ¼ R0 þ dR(l ) (swelling of a Pt _ _ ¼ prl ðR2 dH_ þ 2RHdRÞ into eq. (13) and so crucible) one introduces the correction dM 2 _ _ _ so the _ dH ¼ ð2rl RH 0 =ðrl R  rs r0 ÞÞdR. Since dH is the component in dVC ¼ dV0  dh_  dH, velocity perturbation can be estimated as u2 ¼ dV0 þ dV0R, were dV0R ¼ ð2AHV _ V0 =RÞdR. Since AHV _ < 0 increasing the crucible diameter is equivalent to decreasing the crystal pulling rate and vice versa. In LOM the main perturbations are the melt temperature variations in the meniscus base. The melt temperature depends on the heating power instability and if the heating power is 20 kW to get 1000  C near to the crystallization front, small changes of 20W ((0.1% of heating power) roughly corresponds to 1 C which is a significant quantity. Since it is well known from experiments, that the frequency spectrum of melt temperature fluctuations u1 ¼ dT (t) is very “rich” in the Cz (see work Cariberg [29]), so it is important to find the frequency which is most dangerous in the practical design of feedback control. The presence of pair independent perturbations in the T- and V- channels for some combinations result in an “anomalous” transfer processes at the growing crystal for “normal” densities relations rs/rl > 0. This effect will be presented below. 3.1. Stability analysis of an open state in the Cz method The LOM allows one to make clear the problem of stability in the Cz process since it is necessary to formally answer the question what is the stability of the initial (stationary) state of the growing crystal. We have the freedom to “chose” any arbitrary initial state and carry out a stability analysis of LOM in this open state of the controlled object (Cz). This analysis is based on the determination of the roots of the characteristic eq. (83): DðsÞ ¼ s2  ða11 þ a22 Þs þ a11 a22  a12 a21 ¼ s2  ðspurAÞs þ detA ¼ 0 which always have real-value coefficients. The Rauth-Hurwitz criterion used this algebraic polynomial equation D(s) ¼ 0 and it is not adjusted to dynamic systems, which include a pure time delay [30]. In our case the stability criterion is very simple: spurA < 0

ð95Þ

det A > 0

ð96Þ

It allows one to determine, whether all the roots of a characteristic equation have negative real parts (stability) and if not, how many roots have the positive real parts (instability). It is clearly possible to demonstrate a system stability region with the roots arranged on the complex plane by using the roots locus depending on the technological parameters e the melt

30

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 13. Root locus l1;2 ðDTÞjGs ¼const as function of melt overheating for a LiNbO3 crystal. The first root of the characteristic equation decreases l1 ðdTÞ/  N and does not have any influence on the process stability. But the second (dominant) root quite the contrary increases l2 ðdTÞ/ þ N. In calculations DT ¼ 3 C > 15 C for a predetermined GS ¼ const.

overheating DT and the temperature gradient Gs in crystal at the crystallization front. This stability reserve will be defined by the position of a stable root, nearest to the imaginary axis. A typical Cz root locus, calculated for a LiNbO3 crystal is shown in Fig. 13. The roots of the characteristic equation depend on the material and technology parameters l1,2(DT, GS). The main interest is represented by the dependence of the root values caused by melt overheating l1;2 ðDTÞjGS ¼const . In all cases (for selected typical situations: i) an initial crystallization stage, ii) an intermediate crystallization position at the initial cone and iii) the cylindrical steady state crystal growth) there exists a limit value of the melt overheating for which the system stability region is a maximum (Figs. 13 and 14). An increase over this overheating limit has catastrophic effects by diminishing the stability 10 times from 103 to 104. In practice this means, that the melt temperature variations within

Fig. 14. Every Cz and LEC state has a limited stability region jl1;2 ðDTÞjGs ¼const j ¼ max. This stability reserve is determined by the calculated root l2 ðDTÞjGs ¼const (for the Cz LiNbO3 crystal see Fig. 13 and Table 1).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

31

the limits of several degrees centigrade due to careless operations during seeding or incorrect feedback control in the T- channel can cause stability reduction or even failure, so that the crystal looses contact with the melt (at T w 1500 C (!)). By fixing the melt temperature the stability l1;2 ðGs ÞjDT¼const < 0 improves by an increase of the temperature gradient in the crystal GS (see. Fig. 14). So this analysis, based on LOM shows that the common Cz method has a small stability region, since the typical roots of the characteristic equation D(s) ¼ 0 are small li w 103. And at the same time, the capillarity in the Cz method does not “help” to keep crystallization stability. It is evident in Fig. 15. where the transfer processes h ¼ f(r(t)) are shown, that are caused by a step-wise increase (dT ¼ (2 / 14)  C) and decrease (dT ¼ (2 / 6)  C) of the melt temperature. Stationary cylindrical crystal growth corresponds to the point A(r0, h0) on the curve of the stationary meniscus height h0(r0, 3). In the case where the temperature increment at the meniscus base dT > 0 the crystal can not “find” a new stationary position and loses its contact with the melt where the meniscus height rises to its limiting height Hmax (ri). This height Hmax (ri) is determined from the Jacoby equation discussed below (see Section 8.4). So the dynamic stability of the Cz method is determined not by capillary, but only by the thermal conditions. The physical condition which fulfils eqs. (95) and (96) consists in the presence of any melt overheating at the meniscus base. Really, for the overheated melt the crystallization front tries to return to an initial equilibrium state since below it is hot, and above it is cold. And it is evident that positive thermal perturbations dT > 0 are more dangerous to the crystal having small temperature gradients.

Fig. 15. Transfer process for a Ge crystal in an open state (calculated by use of LOM eqs. (87) and data in Table 2) for a step-wise increase and decrease of temperature DT at the meniscus base. The perturbations are from dT ¼ 6  C to dT ¼ 14  C. An initial stationary crystal growth with r0 ¼ 10 mm corresponds to point A (r0, h0, e). The stationary meniscus heights h(r0, e) are calculated using the Tsivinskii approximation (see part 2.1). The maximum meniscus heights H max(r) were calculated as in [19]. Here we used the corresponding approximation equation where Hmax ¼ Ur, where U(r/a) ¼ 0.4859(r/a) þ 0.2474 for (r/a > 0.2), (see Section 8.4). The crystal separation from the melt would occur at a melt heating of dT ¼ 12 C.

32

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 16. Typical amplitude and phase-frequency characteristics (Bode diagrams) for T- and V- independent perturbation channels. The calculations used LOM eqs (96) and LiNbO3 parameters for an ordinary technological Cz process (see Table 1). (a) stage of crystal pulling r(t) ¼ (0.2  0.5  0.95)r0 (b) T-channel and (c) V-channel perturbations.

3.2. Phase frequency characteristics for Cz in the open state The LOM eq. (87) and standard Matlab library sys ¼ ss(a,b,c,d ) and bode(sys) allow one to calculate the phase-frequency characteristics for the T- and V- independent perturbation channels, the so-called Bode diagrams. In Fig. 16 are shown the calculated results that are typical for the Cz method. It is fundamentally important, that the temperature perturbations dT( f ) at frequencies f > 103 Hz are efficiently suppressed. In other words the growing crystal is represented as a low-frequency filter through which only ultra low frequency perturbations f > 103 Hz pass in the T- channel with a frequency cut-off f z 104 Hz which weakly depends on the crystal diameter. Filtering such ultra low frequency perturbations in the weight measurement signal using a convenient approach represents considerable complexity [31] and the Section 5.1 is specially devoted to a solution of this problem using digital Kalman filters. The V-channel of pulling rate perturbations is practically “transparent” for all noise frequencies. This is obvious in the physical sense, since the velocity perturbations directly

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

33

deform the melt meniscus and such a deformation can include all frequencies. At a frequency f w 103 Hz the Cz process shows oscillatory characteristics, which are more evident, the greater the values of the temperature gradient in the growing crystal. Similar behaviour occurs for Cz growth of different crystal materials. 3.3. Czochralski input response To complete the open-loop analysis of the dynamic Cz system, it is necessary to consider the controlled object output response (transfer processes) to different controllable inputs. The main external input perturbations can be described by a vector U(t) ¼ [dT(t), dV0(t)]T ¼ [u1(t), u2(t)]T with components dT(t), the melt temperature, and dV0(t), the crystal pulling rate variations. They will be approximated by ordinary step functions u(t) ¼ 1(t), having different signs and amplitude. Here we will survey the problem of the so called “anomalous” system response, which is widely considered in the literature since the publication of the well-known work of Bardsley and co-authors [32,33]. The transfer process calculations are based on the discrete linear Cz model and the well known calculation technique (see the book Kuo [34]). The crystallization rate variation can be defined from the eq. (87) set of LOM in discrete time (Section 2.6): 8 _ _ < dVc ðkÞ ¼ dV0 ðkÞ  dhðkÞ  dHðkÞ _ ¼ c11 drðkÞ þ c12 dhðkÞ þ d11 dTðkÞ þ d12 dVðkÞ YðkÞ ¼ prl R2 dHðkÞ : _ dhðk þ 1Þ ¼ f12 drðkÞ þ f22 dhðkÞ þ B12 dTðkÞ þ B22 dVðkÞ as: dVc ðkÞ ¼ dV0 ðkÞ  RRdrðkÞ  HHdhðkÞ  TTdTðkÞ  VVdVðkÞ where RR ¼ ðf21  c11 =PKÞ; HH ¼ ðf22  c12 =PKÞ; ðB22  d12 =PKÞ, PK ¼ prl R2 and R is a crucible radius.

TT ¼ ðB21  d11 =PKÞ;

ð97Þ VV ¼

3.3.1. Response to perturbations having the same sign A single-channel perturbation corresponds to the primary situation, when the system input is stimulated by either temperature dT(t) or pulling rate dV0(t) variations only. First consider the dynamic response for materials which have the density relation rs/rL > 1, i.e. the solid phase is heavier compared to the liquid phase (example: LiNbO3). Second, consider the dynamic response for simultaneous vectored input perturbations for the temperature and crystal pulling velocity. The calculated responses are shown in Fig. 17. The curves 2 and 3 represent single channel (scalar) perturbations dT or dV0 and curve 3 represents a transfer process with a vectorial input perturbation dU ¼ [dT, dV0]T. At the positive scalar (single-channel ) perturbation of the melt temperature dT ¼ þ4  C or the pulling rate dV0 ¼ þ12% (from V0 ¼ 8 mm/h) the crystal diameter decreases (curve 2,3 in Fig. 17b) and the meniscus height increases (curve 2,3 in Fig. 17c). As a result of such perturbations the crystallization front “runs away” into a colder zone and at the end of the transfer process the system comes to any new stationary state. But note here, that for positive scalar pulling rate variations dV0 > 0 the total weight signal Y ¼ dm_ þ dm_ is positive at the beginning of the transfer process (curve 2 in Fig. 17f) and is always negative for the scalar positive input temperature dT > 0 (curve 3 in Fig. 17f). Such an increase in the output signal Y (t) > 0 while the crystal diameter decreases, dr(t) < 0 for the material with the density relation

34

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 17. Response of a LiNbO3 crystal to scalar and vectorial perturbations. The external perturbation has the same (positive) sign (dT ¼ þ4 C and dV0 ¼ þ12%). The initial state is shown in figure (a). We use a crucible size Rcr ¼ 100 mm. (d) At dT ¼ þ4  C and dV0 ¼ þ12% the crystal decreases its radius dr(t) < 0, but the measured signal at the same time is positive Y(t) > 0. (h) Response is mirror symmetric for a LN crystal with the same technological parameters but with the negative perturbations dT ¼ 4  C or dV0 ¼ 12%.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

35

Fig. 18. Transfer processes in a LiNbO3 crystal calculated at scalar and vectorial perturbations. dT ¼ 2  C, dV0 ¼ 22% (V0). Crystal initial state is shown in Fig. 17(a). The output weighing signal Y ¼ dm_ þ dm_ decreases (see (c)) while the crystal diameter increased (see (a)).

rs/rL > 1 should not be considered in any way “anomalous” or irregular. It is the natural consequence of positive pulling rate variations DV0 > 0. It is possible to calculate the meniscus mass variation dmðtÞ _ (Fig. 17g) separately from the _ variation of the crystal mass rate dmðtÞ (Fig. 17f). These two aspects include the summary _ measured by a weighing sensor (Fig. 17d). For positive pulling rate weight signal Y ¼ dm_ þ dm, perturbations dV0 ¼ þ12% in this example the crystal mass rate exceeds the meniscus mass

36

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

_ variation dmðtÞ > dmðtÞ, _ so that the contribution of the crystallization rate variation dVc ¼ dV0  dh_  dH_ in the crystal mass variation dm ¼ prsr0(r0dVc þ 2Vc0dr)dt is more important than decreasing the crystal diameter dr(t) (curve 2 in Fig. 17b). It is evident that the situation is mirror symmetric for a crystal with the same technological parameters (see Fig. 17h) but with negative perturbations dT ¼ 4  C or dV0 ¼ 12%. So for the common Cz with weighing control with the scalar output signal Y(t) and vectorial input perturbations dU ¼ [dT, dV0]T, having the same sign plus (þ) or minus (-) crystallization of materials, having a “normal” density relation (i.e rs/rL > 1) it is possible to register an “anomalous” behavior of the output signal (Fig. 17d). For practical feedback control design that means, the feedback control system, optimized for use of weight signal Y(t) can incorrectly respond if a puller generates a vector input perturbation. 3.3.2. Response to vectorial perturbations having a different sign An “unusual” system response Y(t) should be expected for the materials having the “normal” relation of densities rs/rL > 1 at vectorial perturbations dT(t) and dV0(t) with the different signs and absolute values (Fig. 18). Here the response is given by the fundamental principle of superposition in a linear system and is the sum of each system component. As a result, there exists a possibility of incorrect action of a feedback controller if its design is incorrect. So, if a PID controller has a singlechannel T- feedback, the last situation is “unusual”, i.e. the controller does not “await” the presence of velocity perturbations (dV0 ¼ 22%). Its incorrect action in this case consists in decreasing (not increasing) the heating power due to a decreasing output signal Y(t). Such power reduction has inevitable consequences resulting in an additional increase in crystal diameter, i.e. over control. In general, the Cz gives birth to different transfer processes, but only for outputs due to incorrect assessment of real puller perturbations is it possible to consider as “anomalous”. 3.3.3. Response for material, having rl/rs > 1 In the publications of Bardsley et al. [32,33] the particular role of the meniscus in weight control was discussed for the materials having a density relation rl/rs > 1. But the analysis did not consider the heat transfer issues. Now it is possible to calculate the transfer functions for materials with such a density relation, for example for germanium (Ge), gallium arsenide materials (GaAs) and other materials, having such a density relation rl/rs > 1.The meniscus weight variations are especially important in crystal diameter control based on weighing control for the materials with rl/rs > 1. In Fig. 19 are shown the calculated transfer processes for scalar temperature perturbations dT ¼ 6 CO14 C (at dV0 ¼ 0) for germanium (Ge, Table 2). The output signal is always “singular” at the beginning of the transfer process (Fig. 19d) since the main contribution to the total weighing signal Y(t) made the meniscus weight variation dmðtÞ. _ Note, in this calculation some changes of the crystal front shape were not considered. But for such materials special attention should be given to variations of the crystallization rate dVcr(t) (see Fig. 19f) and eq. (98) and its dependence on crystal radius variations as function of its length dr(lcr) (see Fig. 19g). Formally, calculations can show melting of a growing crystal dVcr(t)/Vc0 < 100% for positive melt temperature changes dT(t) > 0 (the length of growing crystal decreases) which is possible in principle. However the calculated variations are so great, that it is necessary to make corrections to the assumptions in our simple calculation of the transfer processes.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

37

Fig. 19. Calculated transfer processes for a Ge crystal for scalar perturbations dT ¼ ð6O14Þ C (at dV0 ¼ 0). Crystal radius in (a) corresponds to the initial state A in (b). In (f, g) is selected the area in which the crystal formally should remelt during the response.

38

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Table 2 Parameters for the Ge crystallization process Technological parameters (well-known) Seed crystal radius Crystal radius Radius of crucible Angle of crystal cone Crystal pulling rate

(m) (m) (m) (grad) (m/sec)

Rs R0 Rcr b V0

2.5*0.001 10*0.001 20*0.001 50 4*(0.001/3600)

Ge material parameters Thermal conductivity solid phase Thermal conductivity liquid phase Latent heat Growth angle Capillary constant Density of liquid phase Density of solid phase

(W/mK) (W/mK) (J/kg) (grad) (m) (kg/m3) (kg/m3)

ls ll L e0 A rl rs

17.4 9 478000 15 4.77/1000 5.51*0.001/(0.01)^3 5.26*0.001/(0.01)^3

DT GS

4 1500

Technological parameters (rough estimate for Ge crystallization process) Melt overheating (K) Temperature gradient in crystal (K/m)

To a the first approximation it is assumed, that the low bound of the crystallization rate is zero, i.e. dVcr(t) ¼ Vc0. In other words, we shall prohibit remelting of the pulled crystal. At the same time, when dVcr(t)/Vc0 < 100%, the crystal stops increasing its length, but continues its movement upwards at a constant crystal pulling rate V0 ¼ const. Fig. 20 shows the calculated transfer processes for the same germanium (Ge) crystal growth conditions and the same single-valued perturbations T ¼ 6 CO14 C, dV0 ¼ 0, but allowing the above restriction dVc  VC0, i.e. the crystal melting is excluded from exceeding the transfer processes (compare with the Fig. 19). Comparison of the results of calculations (Figs. 19 and 20) reveal that in the second case the melt meniscus does not move up to its maximum height Hmax (curve 2), at which detachment occurs for the same temperature perturbation (dT ¼ 12  C). When crystal remelting is prevented, the increase in height is slower than when the crystal is pulled upward h(t) ¼ V0*t from the melt. It is obvious that in this case, the lateral crystal surface (in our calculations on curves in coordinates rcr(lcr)) has sections with sharp edges. The mathematical model enables one to calculate the dynamic response to periodic perturbations. Fig. 21 shows the response for a Ge crystal for dT(t) ¼ dTsin(2kp/N) with three different amplitudes. Similar perturbations were observed in the work of Hurle et al. [35]. The practical importance of having the estimation of the magnitude of the pulling rate or heating power changes (at the final stage of crystal pulling) which cause the crystal to detach from the melt is that this event is not seen visually and is poorly sensed by a weight sensor. The LOM provides a quantitative estimate of these input variations, which guarantee the crystal detachment with minimum external operations and lower thermal stress. This is especially important for the quality of ferroelectric crystals. In Fig. 22 we estimate the detachment conditions for step-wise changes for the LiNbO3 crystal. 4. Analysis of a closed loop state for Cz Analysis of the open loop state shows that the Cz method belongs to the static controlled object class, having small stability regions. This means, that small changes of melt temperature

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

39

Fig. 20. Dynamic response of a Ge crystal for the same conditions as in Fig. 19. The calculations take account of the remelting restriction for the Cz crystal, i.e. always dVc  VC0 (see also the eq. (98)).

or pulling rate causes irreversible changes of crystal diameter and meniscus height. Formally, by use of the final value theorem [30] applied to the LOM transfer functions Wij(s), defined in eq. (78) we can obtain the open state limðsWij ðsÞuðsÞÞjs/N ¼ const. Here u(s) is the Laplace transform of unit step functions, i ¼ r, h and j ¼ T, V. To avoid the steady state errors in Cz growth it is necessary to use feedback, which has integral components in the T- and V-controlling channels. First of all it is necessary to examine the class of parametric (PID) controllers with T- and V-channels, and in particular, there is great value in having the opportunity to use two control channels jointly.

40

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 21. Response of a Ge crystal to the harmonic melt temperature perturbation dT(t) ¼ dTsin(2kp/N). The initial crystal state corresponds to the point A in Fig. 20a. The calculation takes into consideration the remelting restriction of the Ge crystal during Cz growth. When the amplitude of melt temperature oscillation dT(t) > 4C the variation of crystallization velocity occurs as dVcr(t)/Vc0 < 100% so the crystallization rate Vcr ¼ Vc0 þ dVc reduces to zero (i) and the crystal stops increasing in length and moves up as a massive whole body with a pulling rate V0. In this case the crystal surface develops sharp edges (see also photo in Fig. 50b).

4.1. Determination of an optimum structure for feedback control with parametric regulators. The closed-loop scheme for the Cz system is shown in Fig. 23. The feedback control can be represented by transfer functions WfT(s) and WfV(s) in the T- and V-channels.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

41

Fig. 22. Estimation of the conditions of LiNbO3 crystal detachment from the melt using a LOM calculation for different values dT s 0 and/or pulling rate dV0 s 0 variations (parameters in Table 1). (d,g) e The crystal separation from the melt takes place when the meniscus height is hðtÞxHmax . (b,c,e,f) e with a crystal shape before it separates from the melt. In such an experiment it is possible to identify the actual parameters for a Cz and a LEC dynamical system.

The T-channel input consists in a change of melt temperature after some heating power changes. The melt temperature response following momentary step-wise heating power DN(t) changes has considerable inertia. The most typical temperature response is approximated by the dependence DTðtÞ ¼ ku ð1  expðt=tÞÞDNðtÞ, which is typical for the aperiodic first order

Fig. 23. The block-scheme of a controlled object in signal space. The dotted line shows the feedback transfer functions in the T- and V- channels. The temperature response dynamics for the T-channel is given by the transfer function WT. (a) e block-diagram for Cz with scalar perturbation (in T-channel only); (b)- the general block-scheme for Cz with to take into consideration two perturbations and two control e the T- and V-channels in closed state.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

42

system [30] and is represented in the T- channel of LOM by an additional transfer function WT(s): ku DNðsÞ ¼ WT ðsÞDNðsÞ ð98Þ u1 ðsÞ ¼ DTðsÞ ¼ ts þ 1 Here ku ¼ DTðt/NÞ=DN0 is the gain of this first order system (a coefficient, depending on the heating zone design) and t - is the time constant which is different for heating and cooling of the crucible (more accurate, for the cases of increasing and decreasing heating power). The response (98) is typical for induction heating but for indirect crucible heating there exist more complex temperature responses which are described by the more complex function WT(s) and includes several standard first and second order dynamical systems. It is easy to get the transfer functions for the input perturbations or controlling vectors g(t) ¼ [DT,DV]T and the output error signal vector 3 ¼ (31, 32)T. In the following analysis it is enough to use the unit step approximation for the components of vector g(t). By use of Fig. 23 we obtain the equations set: 8 < Y ¼ WYT 31 þ WYV 32 31 ¼ DT  WT WfT Y : 32 ¼ DV  WfV Y

ð99Þ

then 31 ðsÞ ¼

1 þ WfV WYV 1 þ WYV WfT WT g1 ðsÞ  g2 ðsÞ LðWÞ LðWÞ

ð100Þ

32 ðsÞ ¼

1 þ WfV WYT 1 þ WYT WfT WT g1 ðsÞ þ g2 ðsÞ LðWÞ LðWÞ

ð101Þ

YðsÞ ¼

WYT WYV g1 ðsÞ þ g2 ðsÞ LðWÞ LðWÞ

ð102Þ

where gi(s) e is the Laplace transformation of the input vector components, L(W ) ¼ 1 þ WfVWYV þ WTWfTWYT e is the characteristic polynomial, which determines the stability dynamical system in the closed state, WfT(s) and WfV(s) - are the required transfer functions for feedback in the T- and V- controlling channels in Fig. 24. The system state variables X(s) ¼ [r(s), h(s)]Tcan be determined from the equations: XðsÞ ¼ WiT 31 ðsÞ þ WiV 32 ðsÞ; i ¼ r; h

ð103Þ

All transfer functions in eqs. (100)e(103) were determined in Section 2.5 except WfT(s) and WfV(s). First of all it is necessary to determine an optimum feedback structure by use of the limiting theorems and the determination of error coefficients [30]. It is possible to determine suitable controller structures Wfj(s), j ¼ T, V by establishing the property limit limðs3i ðsÞÞjs/0 ¼ 0. This task was solved by selecting in eqs. (100)e(103) suitable transfer functions Wfj(s) from the class of three-parameter PID regulation rules and the pair of available inputs g(t) ¼ [DT,DV]T. It is known that for a classical PID controller, the parameters that allow one to achieve the best set point performance generally differ from the parameters that allow one to achieve the best reduction of external disturbances. To determine these parameters simultaneously, it is necessary to use PID controllers with two degrees of freedom (see articles Denisenko [36e39]).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

43

Low-frequency temperature perturbations always exists in Cz practice due to slow melt level reductions in the crucible. Therefore an optimized controller should guarantee at least a definded state in the T-controlling channel. After substitution of the transfer functions: WYT ðsÞ ¼ sWHT ðsÞ ¼ AðsÞ=DðsÞ WYV ðsÞ ¼ sWHV ðsÞ ¼ BðsÞ=DðsÞ 1=t WT ðsÞ ¼ ku =ðs þ aÞ; a ¼ Wfj ðsÞ ¼ kDj s2 þ kPj s þ kIj s; j ¼ T; V

ð104Þ

into eqs. (100)e(103) and evaluating in the limit of limðs3i ðsÞÞjs/0 one defines the feedback controller structures in the presence of single- or two-channel perturbations. In general, when the pair perturbations act in a system it is necessary to determine the static errors 3i(0) ¼ Eijgj(0), (i ¼ 1,2; j ¼ 1,2) for different designs of the two-channel feedback controllers. Here Eij e is the (2 2) matrix of coefficients, which are the free items of the transfer functions in eqs. (100)e(103). Table 3 shows the results of such an analysis. In article [40] it is shown that the best Cz feedback dynamics in the case of one-channel perturbation g1 ¼ dT(s) is achieved using two-channel control (the full PID in T-channel and only P in V-channel) since Eij ¼ 0 (see Table 3). However, with this structure it is impossible deliver zero static error if the second perturbations g2 ¼ dV(s) s 0 since: lims31 ðsÞgi ð0Þ=sjs/0 ¼ ðB0 =A0 Þg1 ð0Þ lims32 ðsÞgi ð0Þ=s ¼ g2 ð0Þ where B0 ¼ Arr Ahh AHV  Arh Ahr AHV þ Arh AHr and A0 ¼ Arh AhT AHr . In practice it is possible to make these state errors small if one keeps DVjt/N ¼ g2 js¼0 /0, so that the melt temperature variations are the main perturbations. In this case, only the P- or D- components are needed in the V-control channel. In the presence of I-integral action in the V-control channel then the Table 3 Values of the static coefficients Eij for the error matrix for various structures of feedback control at standard input vector perturbations gi(s), (i ¼ 1,2; j ¼ 1,2).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

44

pulling rate result in a non-zero error equal to lims31 ðsÞg1 ð0Þ=sjs/0 ¼ ðkIV B0 a=A0 Þg1 ð0Þ. This static error kIV ¼ 0 if the I-component in the V-controlling channel is absent. These conclusions may be extended to the case of nonstationary processes for the growth of the initial (or final) crystal cones, because the structure of the transfer functions do not change. 4.2. Initial and final crystal cone shapes for Cz using weight control The task consists in forming an initial cone with a predetermined shape (this is especially b which important for the LEC method). It is necessary to determine the nonlinear program YðtÞ is monitoring the feedback controller, so that at every regular controlling event it compensates any input step dY(t) similar to any input perturbations g(t). The weighing control program task can represented as r ¼ r(z), b ¼ b(r) or b ¼ b(z) as a function of time, were r, b-are the current values of the crystal radius and slope to its lateral surface at the three-phase line. The total mass for the common basic Cz process measured by a weighing sensor at every time t is: Z t 2 _ ð105Þ M ¼ prl Rcr HðtÞdt 0

_ Where HðtÞ is the rate of melt level decrease in the crucible with radius Rcr ¼ const. Otherwise: Z t ð106Þ M ¼ prs r 2 ðtÞVc ðtÞdt þ rl Vm  ðrs  rl ÞUs þ m0 0

_ b; tÞ  HðtÞ _ Here Vc ¼ V0  hðr; is the crystallization rate, V0¼const is the crystal pulling 2 2 rate, Vm ðtÞ ¼ prl r ðtÞh þ a prl rðtÞcosð30 þ bÞ- the meniscus volume, a ¼ ð2sSL =rl gÞ1=2 is the capillary constant, Us(r,t)-the volume of the crystallization front deflections (plus corresponds to a convex and minus to a concave shape of the crystallization front, see Fig. 6) and m0¼const is the mass of seeding crystal and its holder. For the determination of the crystallization rate Vc(t) it is necessary to specify an appropriate function b(r) and to calculate the derivatives:   _ ¼ M 0 þ M 0 þ M 0 b0 r_ M t r b r Mt0 ¼ prs r 2 Vc Mr0 ¼ 2prl hr þ prl r 2 h0r þ prl a2 cosð30 þ bÞ  ðrs  rl ÞU0sr Mb0 ¼ prl r 2 h0b  prl a2 rsinð30 þ bÞ   H_ ¼ V0  Vc þ h0 þ h0 b0 r_ r

ð107Þ

b r

Here derivatives h0r , h0b were derived using eqs. (5) and (6) in Section 2.1. Taking into account the relation r_ ¼ Vc tgðbÞ and the equation set (107) it is possible to determine the value of the variable crystallization rate for every time event of the program task. For the simplest approach of a flat crystallization front Us¼0 or more exactly, when the change in the rate of the front _ s /0, we get [40]: volume is a negligible quantity during initial crystallization of the cone U Vc ¼

pr R2 V  l cr 0   i prl R2cr  prs r 2 þ prl R2cr h0r þ h0b b0r  Mr0 þ Mb0 b0r tgðbÞ h

ð108Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

45

Fig. 24. (a) A typical nonlinear program task for the initial crystal cone (photo LiNbO3 crystal, left) using the predetermined function angle changes b(r). (b) a similar approach for the final cone in the Cz method.

The last assumption is quite appropriate for most materials pulled from the melt by the basic Cz process. For r/r0 and b/0 eq. (108) turns into the well known Vc0 ¼ rl R2cr V0 =ðrl R2cr  rs r02 Þ. The behavior of the curvature of the crystallization front depends on the crystal material, thermal conditions in the growth zone and the melt hydrodynamics in the crucible. Calculating the value Us(r,t) in controlled (real) time is difficult today or even impossible due to the complexity of this mathematical task. But it is possible to do quite efficient approximations, based on a posteriori processing of the experimental results. The simplest approach may be based on the use of a sphere segment, having a variable curvature radius Rfr(t), as the shape of the crystal front. The time dependence Rfr(t) may be successfully estimated experimentally and can be used in the program task calculations with the set up of a feedback controller. In practice growing an oxide crystal, (for example LiNbO3), the curvature is small R1 fr ðtÞ/0, so it is correct to assume that the crystallization front is flat Us¼0. The algorithm for the calculation of the program task in discrete time is: -at i-th step using the predetermined dependence b¼b(r) we find bi ¼ bðri Þ -compute the value of the current crystallization rate Vci ¼ Vc ðri ; bi Þ (see eq. (108)) _ i ; bi ; Vci Þ ¼ _ i ¼ Mðr -determine the program rate of change of the total measured weight M _ prl H i ðri ; bi ; Vci Þ (see eq. (107)) -prepare the following step riþ1 ¼ ri þ T0 Vci tgðbi Þ and calculate the subsidiary integral values _i Miþ1 ¼ Mi þ T0 M liþ1 ¼ li þ T0 Vci Using the calculated current crystal length li it is possible to reconstruct the cone shape in the usual coordinates r¼r(z). Fig. 24 shows the practical results for the initial and final cones based on the use of the dependence bðrÞ ¼ bm  gðr  r1 Þ2 , were bm is the greatest slope angle to the crystal cone profile curve, g and r1 are the parameters, derived from the boundary conditions

46

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

b(0)¼bs and bðrfin Þ ¼ b0 , were rfin is the final cone radius, and bs and b0 are the initial and final slope angles of the final crystal cone. If it is necessary to calculate a predetermined cone shape as function of its length one can use any appropriate dependence b(z). This problem is especially important in the automated Liquid Encapsulated Czochralski method and will be discussed in Section 6. As an example consider a very simple and useful dependency bðzÞ ¼ bm sinðgzÞ for which the condition for the transition from a cone to a cylinder is gzc¼p. Here zc is the required crystal cone length. In order to represent the relation r_ ¼ Vc tgðbÞ as dr=dz ¼ tgðbÞ so: Z zc   tg½ðbm sinðpz=zc ÞÞdz ¼ zc 30bm þ 4b3m 45p ð109Þ Dr ¼ r0  rseed ¼ 0

P 2n This integral is obtained following expansion of the value tgðxÞ ¼ N n¼1 2 2n1 ð2  1Þ=ð2nÞ! Bn x in the Taylor series at jxj < p=2, B1¼1/6, B2¼1/30 [41]. From this equation it is possible to determine the crystal cone length zc ¼ 45pDr=ð30bm þ 4b3m Þ and the parameter g ¼ p=zc. The initially predetermined parameters used were: - crystal radius r0, seed radius rseed and bm- the largest slope angle to the cone profile. Further computations are similar to those represented by eq. (107), taking into account the substitution b0r by b0z. In general, the necessity for quadrature is restricted by the cone shape according to the choice of an appropriate function b(z). Also it is possible to get the predetermined cone shape r¼f(z)or z¼f(r) by a similar method, since b is determined from ðdr=dzÞjri ¼ tgðbi Þ or ðdr=dzÞjzi ¼ ctgðbi Þ. And at last it is possible to calculate the desired cone shape using parametric functions, for example, as r¼r(t), z¼z(t). 4.3. Optimized feedback control tuning (Parametric regulators) Examples of transfer processes for different PID control gains for single- and two-channel PID feedback is given in [40]. In a discrete model the temperature response for a similar step power change can be described as DTðk þ 1Þ ¼ 4DTðkÞ þ kU ð1  4ÞDNðkÞ

ð110Þ

where 4 ¼ expðT0 =tÞ, T0- is the sampling time. So such an approach allows one to use only the “natural” state variables, which have a definite physical sense e the melt temperature, crystal radius andRmeniscus height variations. With standard PID control in continuous time t uðtÞ ¼ kp ½eðtÞ þ T1I 0 eðxÞdx þ TD e_ ðtÞ it is possible to rewrite in recurring form in discrete time uðkÞ ¼ uðk  1Þ þ q0 eðkÞ þ q1 eðk  1Þ þ q2 eðr  2Þ, where q0 ¼ kp ð1 þ TD =T0 Þ, q1 ¼ kp ð1 þ 2TD =T0  T0 =TI Þ and q2 ¼ kp TD =T0 [42]. It is obvious, that the control in the T-channel initially changes the crucible heating power DNðkÞ and then varies the melt temperature according to the response of a standard First-Order system, described by eq. (100). So, the transfer processes in Cz with PID feedback control can be described by the simple first order equations set: Xðk þ 1Þ ¼ AXðkÞ þ BðgðkÞ  UðkÞÞ YðkÞ ¼ CXðkÞ þ DðgðkÞ  UðkÞÞ DNðkÞ ¼ DNðk  1Þ þ q0 YðkÞ þ q1 Yðk  1Þ þ q2 Yðk  2Þ u1 ðk þ 1Þ ¼ 4U1 ðkÞ þ kU ð1  4ÞDNðkÞ u2 ðkÞ ¼ kpV YðkÞ

ð111Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

47

In eq. (111) we use different PID-law coefficients, an asymmetrical coefficient of melt temperature response inertia t and some combinations of scalar or vectorial perturbations g(k). For the parametrical controller the general controlling quality criterion may be represented as the quadratic-based functional: Z T  T  1 3i Q3i þ UT RU dt ¼ min ð112Þ J ¼ lim T T/N

0

where Q, R eare the (2 2) matrix of weight coefficients, which determine the restrictions on the errors in the T- and V- channels and of the input controlling power. This criterion is correct in general, because the constancy of diameter of the growing crystal is a necessary requirement for the constancy of the growth condition at the crystallization front which determines the crystal homogeneity. But a constant diameter on averageRover a long time provides a constant crystalT lization rate and also with on average Lim T 1 0 Vc ðtÞdt ¼ Vc0 , at the same time there are T/N _ The absolute values of the deviations possible considerable variations dVc ðtÞ ¼ dV0  dh_  dH. may by small, but their derivatives are not small. In other words, the strict achievement of the constancy of crystal diameter can increase the oscillations of the transfer processes and the crystal

Fig. 25. Calculated transfer processes for two channel control (PID)T þ (P)V with an optimum set up of feedback control coefficients. Here the control is produced by every regulator sampling event and does not take into account the noise in the measured signal. The initial LiNbO3 crystal state e cylinder (r0 ¼ 40 mm) at a pulling rate V0 ¼ 3 mm/h with initial perturbations x1 ¼ dr(0)¼2 mm, x2 ¼ dh(0)¼1 mm and x3 ¼ dT(0) ¼ 3 C. Here u1-is the input in the T-feedback control channel, u2-the input in the V-channel and Vc-is the restored crystallization rate.

48

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 26. Photo of a typical series of LiNbO3 crystals which were obtained by automated Cz using an optimum set up of coefficients for the PID feedback controllers without any limitations for the output controlling power (i.e., matrix R ¼ 0 in eq. (112)). The growth bands are visible on the crystal surfaces in this case.

rate variations accordingly. Any rational compromise consists in an informal choice of the coefficients Q and R for the integral criterion value (112), based on a comparison of the crystal quality and its geometry, obtained with different coefficients for the PID laws. It has been shown above that the feedback structure PIDTþPv of feedback control is an optimum in the case where the temperature perturbations in the T-channel are the principal ones (see Table 3). The determination of the optimum kP, TI, TD and kPV coefficients consist in a calculation of a global minimum of the function (112) in four-dimensional space, based on the Powell method [43]. This approach is shown in Fig. 25 with the calculated transfer function for the case of a LiNbO3 crystal, when there are no limitations in the output controlling power (i.e., matrix R¼0 in eq. (112). So, with an optimum feedback control set up in accordance with the quadratic value criterion the closed loop response becomes oscillatory. As a consequence crystals using such “hard” feedback control develop banding segregation. Sometimes it is immediately visible, for example in the LiNbO3 crystals (Fig. 26). Correspondingly for any “intermediate” feedback tuning control, when the coefficient kp in the T- channel was reduced the transfer function oscillations and the melt temperature amplitude also reduced (see Fig. 27). Now take into account the presence of noise in the output signal, which prevents the controller working at every sampling event. If one allows the introduction of feedback control only after n intervals of sampling time T0, in other words it is necessary for the time t¼nT0 for noise filtering, after which the weighing sensor signal dY(t) is free from noise before the next controlling event, so the transfer processes changes compared to the “noiseless” cases shown in Fig. 28. The calculated transfer processes look similar to the previous case (Fig. 27), but all PID coefficients of feedback in the T-channel were considerably reduced as a result it is possible to get good results both in crystal diameter control and internal homogeneity (Fig. 29). The problem of noise filtering will be discussed in Section 5.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

49

Fig. 27. Calculated transfer processes for two channel control (PIDT þ PV) with reduced coefficient kp in the T- channel compare with the “optimum” feedback set up (see Fig. 26). The initial state and the perturbations in LiNbO3 crystal are the same.

4.4. Feedback control in the Cz method using the state variable observer concept Parametric PID controllers do not in general require mathematical modelling of a controlled object to work successfully it is sufficient to use only the signal deviation from the predetermined program values. The resulting simplicity of PID controllers are their main advantage and this allows their wide usage for control in automated crystal growth from the melt. But a (adequate) mathematical model allows one to get much more information about the real process, it enables one to estimate the system state variables and to use modern state controllers for feedback. The potential for the use of state controllers in Cz and LEC methods is discussed below. 4.4.1. Observability, controllability and state controllers in the Cz method With weighing output control (in detail see Section 7) the state variables dr(t) and dh(t) are not measured directly, but can be restored by use of the so called state observer (see the book Kvakernaak, Sivan [44]). This approach requires a preliminary analysis of the problem of the complete operability, controllability and observability of the Cz method, taking into account in general, the presence of stochastic noise in the input and output signals. Writing down the inertia-free LOM of basic Cz in discrete time as the equations set:      a11 a12 x1 ðkÞ b11 b12 u1 ðkÞ x1 ðk þ 1Þ ¼ þ ð113Þ a21 a22 x2 ðkÞ b21 b22 u2 ðkÞ x2 ðk þ 1Þ

50

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 28. Calculated transfer processes for two channel control (PIDT þ PV), which include 6 minute pauses used for noise filtering. The transfer process includes 25 control cycles. The coefficient kp in the T- channel reduced from kp ¼ 1500 to kp ¼ 0.5 compared to the set up of feedback control in Fig. 27. The initial state of the LiNbO3 crystal and the initial perturbations are the same.

Fig. 29. A Photo of a LiNbO3 crystal which was obtained by Cz with PID feedback control and filtering of the weighing measurement noise. Duration of the time interval between crystal diameter control events was 6 minutes. A limit in the output controlling power was taken into consideration (R s 0 in eq. (112)).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

with the output equation (see Section 2.6, eq. (87)):   x1 ðkÞ u1 ðkÞ þ ½ d11 d12  YðkÞ ¼ ½ c11 c12  x2 ðkÞ u2 ðkÞ

51

ð114Þ

Here x1(k)¼dr(k) and x2(k)¼dh(k) are the state variables and variations u1(k)¼dT(k), u2(k)¼dV0(k)- are the components of the input vector of the Cz system. If the initial state of the system is represented by steady cylindrical crystal growth (b(t)¼0) and correspondingly b11¼b12¼0), so the structure of this simplest model (113) and (114) for the open loop Cz method looks like: Since operability and controllability are the structural properties of a dynamical system, represented by its mathematical model [27], so it follows directly from this decomposition scheme that the output signal of the weighing sensor Y(k) “holds” the information about both state variables e the meniscus variations x2(k)¼dh(k) and the crystal diameter variations x1(k)¼ dr(k). Formally this is confirmed by the Kalman criterion [34]: e the rank of the observability matrix ½CT jAT CT  ¼ 2 is equal to the dynamic system order. Also directly from Fig. 30 it is obvious, that the basic Cz method is completely controllable since both inputs (in the T- and V- channels) have an effect upon the meniscus x2(k) and the crystal subsystem x1(k). And this also formally conforms with the rank of the controllability matrix ½BjAB ¼ 2 which is also full. Testing the LOM obtained by using standard Matlab programs Co¼ctrb(A,B) and rank(Co). It is easy to see that complete controllability of the Cz method is possible even with one-channel input. However it is necessary to remember, that the verification of the system observability and controllability does not provide an answer to what control quality is achievable. So, for example, it is evident that it is not possible to compensate for every melt temperature perturbation dT using single V- channel feedback control of the Cz process. And since the pulling rate perturbations in practical Cz are either very small or completely absent T- single channel feedback control can be successfully used, based only on melt temperature changes. Following on from the material discussed above suppose that the more realistic mathematical model must additionally include the melt temperature inertia model as a response to the step wise heating power changes. So an improved model for the controlled object (Cz) can be rewritten as an expanded one: 8 < x1 ðk þ 1Þ ¼ a11 x1 ðkÞ þ a12 x2 ðkÞ þ b11 DTðkÞ þ b12 DVðkÞ x ðk þ 1Þ ¼ a21 x1 ðkÞ þ a22 x2 ðkÞ þ b21 DTðkÞ þ b22 DVðkÞ : 2 DTðk þ 1Þ ¼ 4DTðkÞ þ ku ð1  4ÞDNðkÞ

ð115Þ

Fig. 30. Representation of the basic Cz method with the “inertia-free” Low Order Model in state variables space.

52

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 31. Representation of the basic Cz method using the Low Order Model in state variable space, which takes into account a typical melt temperature response DTðnT0 Þ ¼ ku ð1  enT0 =t ÞDN to the step wise heating power DN changes.

Here DTðkÞ ¼ dT þ dT; DVðkÞ ¼ dV0 þ dV0 are the control inputs and perturbations in the T- and V- control channels,4 ¼ expðT0 =tÞ, T0 is the sampling time, ku DN is the amplitude of step like heating power changes. And the output equation may be rewritten as: YðkÞ ¼ c11 x1 ðkÞ þ c12 x2 ðkÞ þ d11 DTðkÞ þ d12 DVðkÞ Or in matrix notation: 3 2 3 2 3 32 2 a11 a12 b11 4 b11 ku ð1  4Þ b12  x1 ðkÞ x1 ðk þ 1Þ 4 x2 ðk þ 1Þ 5 ¼ 4 a21 a22 b21 4 54 x2 ðkÞ 5 þ 4 b21 ku ð1  4Þ b22 5 u1 ðkÞ u2 ðkÞ 0 x3 ðk þ 1Þ x3 ðkÞ ku ð1  4Þ 0 0 4 3 2  x1 ðkÞ u1 ðkÞ 5 4 YðkÞ ¼ ½ c11 c12 d11 4  x2 ðkÞ þ ½ d11 ku ð1  4Þ d12  u2 ðkÞ x3 ðkÞ

ð116Þ

ð117Þ

ð118Þ

Here the third equation in (118) describes the discrete form of a typical (and the simplest) melt temperature response DTðnT0 Þ ¼ ku ð1  enT0 =t ÞDN. The full observability and controllability also conform formally since the rank of the observability matrix ½CT jAT CT jðAT Þ2 CT  ¼ 3 is equal to the order of this expanded system and the order of the controllability matrix ½BjABjA2 B ¼ 3 is also full. Similarly, the complete observability and controllability may be established from the decomposition scheme in Fig. 31.

Fig. 32. (a)-Controlled object (Cz) can be used both with the PID controller in feedback, and the state controller, based on a Luenberger observer. In such a configuration the state observer can only restore the system state during the work of the digital PID controller in feedback mode. In this case the matrix coefficients K ¼ 0. Detail in the text; (b) e representation of the Cz process in discrete closed loop with weight output control in state space.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

53

Note here, that the model (117) and (118) makes sense if the input (feedback) control uses both the inputs dN and dV. In this case dN(k)¼0 and dTis the inertia free melt temperature perturbations at the meniscus base. But if one only uses the V-channel it is sufficient to utilize the inertia free Cz LOM process. The full observability and controllability conditions and the authentic information about the input control u1(t) ¼ dN(t) and pulling rate u2(t) ¼ dV0(t) allows one to use the full Luenberger observer for the reconstruction of the Cz state variables. With weighing control the state variables in the Cz system are not measured directly but can be reconstructed and then used in the optimum state controller. A scheme for such a combination of state controller and observer is shown in Fig. 32. It is possible also to use a PID controller in closed loop Cz and reconstruct simultaneously with the state observer for the state system variables. This possibility will be discussed below. In this approach the state observer uses the deviations from the system steady state as the initial state of the dynamic system (LOM), which should be reconstructed using the Luenberger observer concept [42,44]. For system observability it is necessary to establish the following conditions: i) the mathematical model (117) and (118) must adequately describe the real object, ii) the object parameters are known, such as the time constant t of the temperature response for the heating power changes, the melt overheating at the meniscus base and the initial temperature gradient GS jr¼r0 in the growing crystal near its crystallization front, iii) correct information is available about the input control for the heating power u1(t) ¼ dN(t)and pulling rate u2(t) ¼ dV0(t). Of course, the output signal from the weighing sensor must be free from measurement noise. The equivalent observer can be represented by the equation set:  b b þ 1Þ ¼ A XðkÞ þ B UðkÞ þ HDe Xðk ð119Þ b YðkÞ ¼ C XðkÞ þ D UðkÞ b Here * marks the expanded matrices in eqs. (117) and (118), De ¼ YðkÞ  YðkÞ and b b XðkÞ; YðkÞ are the calculated values for the state variables and the observed output values of the weight sensor. The accuracy and the convergence rate of the observer is given by b b b ¼ 0, where DXðkÞ ¼ XðkÞ  XðkÞ are determined by choice of the roots lim ðA  HCÞDXðkÞ k/N  of the characteristic equation detðzI  A þ HCÞ ¼ 0: 3 2 z  a11 þ H1 c11 a12 þ H1 c12 b11 4 þ H1 d11 4 ð120Þ det4 a21 þ H1 c11 z  a22 þ H2 c12 b21 4 þ H2 d11 4 5 ¼ 0 H2 c12 z  4 þ H3 d11 4 H3 c11 As an example it is possible to use the criterion of quick convergence. In accordance with this criterion the roots li, i¼1,3 in eq. (120) have to be located at the point li¼0 [42] and can be obtained from the algebraic equations set: pij Hi ¼ li p11 ¼ c11 p12 ¼ c12 p ¼ d11 4 Where 13 p31 ¼ a22 c11  a21 c12 p32 ¼ a11 c12  a12 c11 p33 ¼ a11 c22 d11  a12 c21 d11  a22 b11 c11

ð121Þ p21 ¼ c12 a21  c11 a22 þ c11 4 p22 ¼ c11 a12  c12 a11 þ c12 4 p23 ¼ c11 b11 4 þ c12 b21 4  d11 ða11 þ a22 Þ4 l1 ¼ a11 þ a22 þ 4 l2 ¼ a21 a12  a11 a22  ða11 þ a22 Þ4 l3 ¼ a11 a22  a12 a21

54

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

The noise in the output signal considerably restricts the reconstruction time. Since the order of our system k¼3, so theoretically the full reconstruction time is possible after three measurement steps. In practice it is more convenient to determine the root positions inside the unit circumference on the complex plane. This guarantees the convergence of the observability process and decreases the effect of noise, which is unique for every puller. The coefficients of the matrix H, corresponding to the predetermined position of the stable real roots li, i ¼ 1,3 can be obtained from the algebraic equations set: 8 < p1j H1j ¼ l1 þ l1 þ l2 þ l3 p2j H2j ¼ l2  l1 l2  l1 l3  l2 l3 ð122Þ : p3j H3j ¼ l3 þ l1 l2 l3 Here the well known Viet theorem about correlations the of algebraic equation and its coefficients is used. It is possible that this matrix coefficient H, determined from (121) or (122) can be considered as an estimation of the upper boundary values for the all possible observer coefficients. A more general approach for the determination of the observability matrix consists in the use of an integral P criterion. In discrete time it may be represented as N1 T X ðkÞQXðkÞ þ UT ðk  1ÞRUðk  1Þ, were Q and R are the weight J ¼ XT ðNÞQXðNÞ þ k¼0 matrix. As it is well known, the minimization of this quadratic criterion during the determination of an optimal feedback state control UðkÞ ¼ KXðkÞ consists in the solution of the recurring matrix equation set:  1 Knj ¼ R þ BT PNjþ1 B BT PNjþ1 A  ð123Þ Pnj ¼ Q þ AT PNjþ1 A  KT R þ BT PNjþ1 B KNj Here K is the matrix of the final resulting values of such a regulator, PN is the initial matrix and K ¼ lim Knj . N/N By use of the well known duality property of a linear system it is possible after the replacement A/AT , B/CT and K/HT to construct the dual object which may solve the task of the optimal observer. From the solution of the recurring equation set:  1 T HT ¼ R þ CPNjþ1 CT CPNjþ1  A  ð124Þ T Pnj ¼ Q þ APNjþ1 A  HNj R þ CPNjþ1 CT HT Similarly here the matrix PN is the initial matrix and H ¼ lim Hnj . N/N So the dynamic of a Cz system in closed state, based on the LOM eqs. (117) and (118) and the Luenberger observer with the state regulator in feedback (Fig. 32) is determined by the complete set of roots of the characteristic equation [34]: detðzI  A þ BKÞdetðzI  A þ HCÞ ¼ 0

ð125Þ

Correspondingly the simplest state regulator may be represented as the proportional regulator: x1 ðkÞ  ki2b x2  ki3b x3 i ¼ 1; 2 Ui ¼ ki1b

ð126Þ

In accordance with the well-known division theorem [34,42] it is possible to implement the feedback and observer adjustments separately and independently using the solutions of the Ricatti algebraic equations set (123) and (124). The choice of the weight matrix Q and R in an

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

55

Fig. 33. (a) The calculated transfer processes for LiNbO3 growth by Cz in closed loop after step wise (single-channel) b perturbations dT ¼ 2C, dV ¼ 0 and optimum feedback state control UðkÞ ¼ K XðkÞ with 150 cycles of control time T0 ¼ 1 min. Curve 1 ethe system reaction and 2- observer response; (b)- the transfer process with the same initial conditions, but increasing the controlling time interval to T0 ¼ 6 min; for noise filtering the Cz LOM eqs. (117) and (118) were used.

integral criterion J is required, in general, the noise frequency and its level in a real puller is determined experimentally. But in a numerical experiment it is possible to estimate the matrix values of K ¼ lim KNj and H ¼ lim HNj from eqs. (121) and (122). Below in Fig. 33 the N/N N/N transfer processes with some step wise perturbations (x1(0)¼2mm, x2(0)¼1mm, x3(0)¼2C ) are presented emulating the work of real Cz crystal growth, using the state regulator in feedback and the full Luenberger observer reconstruction of the system state variables: A typical situation for a Cz system consists in a melt temperature perturbation (here dT¼ C) and the absence of pulling rate variations. The equivalent observer (119) successfully restored all of the three state variables e the crystal radiusx1 ðtÞ ¼ drðtÞ, meniscus height x2 ðtÞ ¼ dhðtÞ and the melt temperature variation x3 ðtÞ ¼ dTðtÞ at the meniscus base. The state b (eq. (129)) also successfully returned the dynamic system into an equicontroller U ¼ KX librium state. Fig. 33b shows similar control obtained increasing the control interval from 1 min to 6 min, which was needed for noise filtering of the the sensor output (Section 5). It is possible to use the observer concept (119) and the expanded Cz LOM eqs. (117) and (118) for reconstruction of the state variables in feedback using a parametric PID controller (Fig. 32). Really, for the observer it is “unimportant” which type of controller (PID or state controller) is used to change the heating power and the pulling rate at regular intervals since the input UðtÞ ¼ ½dN; dV0 T is measured precisely in the real closed loop system. Fig. 34a shows such a reconstruction of the state vector X(t) with step-wise power change u1 ¼ ku DN ¼ 2:4 C, which can be produced by any type of controller. Fig. 34b shows the supervision of the state vector X(t) by an expanded observer (119) jointly with two-channel (PIDTPv) feedback

56

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 34. (a)-Calculated record of the reconstruction effect of the statevariables xi(t) and crystallization rate Vc(t). This emulation of the Cz growth of LiNbO3 involves an input perturbation dT ¼ 2C (dV ¼ 0). The calculations involve the use of the expanded observer eq. (119); (b)- results for the same object controlled by a two-channel parametrical regulator (PIDT þ PV).

control. This approach allows one to get information about the crystal radius, meniscus height and the melt temperature variations (at the meniscus base), i.e. to reveal the characteristics of the Cz processes in the open and close states (Fig. 34). 4.4.2. A controlled system with random input fluctuations and noise in the output Let us take into account the presence of stochastic perturbations of white noise of type uY(k)in a measured weighing sensor signal Y(k) and self generated system noise uY(k)which exists in the real growth system at every sampling time kT0. The average of the noise distribution for uY(k) and uY(k) is evidently equal to zero and may be emulated in a mathematical experiment using a random number generator. The measurement noise is caused first of all by the imperfection of the puller’s mechanical drive for crystal rotation and its connection to the weighing sensor. Correspondingly temperature fluctuations are caused by the forced convection in the melt. Note, that a white noise approach not accidental but based on a comparison of the noise frequency behavior and the systems own Cz frequencies, which are located in the ultralow range wð103  104 ÞHc (see Section 3.2) and differ by more than 2e3 orders of magnitude compared with the typical highfrequency perturbations. In general the stochastic model of the Cz method can be written as:  Xðk þ 1Þ ¼ AXðkÞ þ BUðkÞ þ UT ðkÞ ð127Þ YðkÞ ¼ CXðkÞ þ DUðkÞ þ UY ðkÞ Here X3 ðkÞ ¼ DTðkÞ, the joint stochastic process ½UT ðkÞ; UY ðkÞ is derived from UT(k) e the melt temperature noise at the meniscus base and UY(k) e the weighing measurement noise is 0 V considered as white noise, having a permanent density V ¼ ½ 1  for which the matrices V1 and 0 V2 V2 have only the components EfuT ðkÞuTT ðkÞg and EfuY ðkÞuTY ðkÞg. These noises are uncorrelated

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

57

amongst themselves and also with the initial state X3(0), having the definite dispersion Ef½X3 ðkÞ  X30 ðkÞX3 ðkÞ  X30 ðkÞT g¼ S0 . The values S0,V1 and V2 can be estimated after the statistical treatment of the experimental data. Again, for the system (127) we can create an observer: b b þ 1Þ ¼ A XðkÞ þ B UðkÞ þ HðkÞ½YðkÞ  C XðkÞD UðkÞ ð128Þ Xðk and can calculate the matrix H, which minimizes the measurement error. As is well known, the solution of the optimum measurement task in the presence of stochastic noises consists in the calculation of the Ricatti matrix equation [42,44]:  A SðkÞ þ SðkÞA þ V1  SðkÞC V1 2 C SðkÞ ¼ 0 The solution S(k) of this equation enables one to get the observer matrix H:

HðkÞ ¼ SðkÞCT V21

ð129Þ ð130Þ

b Similarly, the restored state variables XðkÞ are used in state feedback control b UðkÞ ¼ KXðkÞ. Here the state controller matrix K ¼ lim Kk can be defined as the final k/N resulting values of the equation: KðkÞ ¼ R1 BT PðkÞ

ð131Þ

Fig. 35. (a)- Emulation of the work of an optimum state regulator eq. (131) in the presence of white noise in the sensor output dY ¼ 50%with stochastic characteristics from the experiment (see Section 5) and self generated random melt temperature perturbations of the the Cz process of amplitude dT ¼ DX3 w100%.

58

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

The matrix P(k)- is the recurring solution of the Ricatti matrix equation: A PðkÞ þ PðkÞA þ Q  PðkÞB V2 B PðkÞ ¼ 0

ð132Þ

With the initial condition P(0)¼S0 and matrices Q, R-which are the weighing matrices of integral controlling quality criterion J eq. (112). Fig. 35 shows the calculated transfer processes during the work of such a state regulator in the presence of stochastic noises in the output signal Y simultaneously with the input melt temperature fluctuations. Though there are very considerable noise amplitudes (in 50% for Y and 100% for X3) the observer provides very good control quality. Similar good reconstruction results with the system state variables Xi (i¼1,3) were obtained for simulation work with PID control with Cz feedback using the same noise parameters. The state controller allows one to change the time interval for filtering in different pullers without changing the feedback set up. This is an important advantage for such regulators compared with parametric PID controllers. 5. Problems and solutions for filtering the measurement noise in automated Cz For Cz and LEC using weight control the information about current crystal radius and meniscus height at the crystallization front is derived from the signal deviations of the current _ measured weight rate from the program values Y ¼ DPðtÞ. The differentiation operation increases the noise which always exists in the output signal. So practically the filtering problem is of key importance. With the usual PID feedback control the high noise level increases the oscillation of the transfer processes and can suppress recovery of the state variables by the observer. The capability of analogue (hardware) filtering is restricted down to frequencies of 0.1e0.3 Hz but in Cz and LEC the main noise filtering frequencies are situated in the ultralow

Fig. 36. (a) - The output weight signal for a pulled and rotating Cz LiNbO3 crystal; (b)- the correlation function of this signal; (c) - its spectral density. The principal noise mode is f0 ¼ 0.203Hc. The sampling frequency exceeds by a factor of two the highest signal frequency making this experiment unique [45].

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

59

region (0.001e0.003) Hz, since this is the region where the crystal’s own crystal growing frequencies originate (Section 3.2). Usually the static sensitivity of the weight sensor is better than 10 mg. But crystal rotation can make this value more than 150 times worse, and this value becomes even larger with increasing crystal mass. In Cz the crystal rotation is the main cause of measuring errors and there is always an harmonic component in the output weight signal. A typical sampling output signal is shown in Fig. 36a. The presence of the main harmonics are easy to detect from the correlation function R(t) of this signal (Fig. 36b). It is necessary to filter this non-stationary signal in order that at the timing of the next current input feedback control event only uses the information about its ultralow changes. One obvious approach is to suppress disturbances of this type. Fig. 36a demonstrates a simple moving average algorithm. It is known, that the efficiency of averaging is determined by a sampling length and a dispersion value of a stationary ergodic stochastic process. There exist strictly determined criteria for the precise values for such process statistical events using finite sampling (see book Bendat, Pirsol [46]). The averaging efficiency depends on a noise amplitude and if it is increased, the sampling length nT0 and the time delay of an average signal td ¼ nT0/2 are increased. The second disadvantage of this simple filtering consists in an oscillation of the average value around the real values due to random shifting of a signal phase. The use of a recurrent version of a moving average algorithm is a little more effective and can be represented by the formula byðk þ 1Þ ¼ byðk  1Þ þ 1=k½xðkÞ  byðk  1Þ, were x(k)-is a discrete reading of a sensor signal and byðkÞ -is a restored average value. Similarly, in this case it is also possible for the oscillation of the average value, which is situated in an information band to cause forced oscillations of the control input and crystal diameter. The attempt at noise filtering in the frequency range near to 0.001 Hz by use of this moving average method causes an general the digital filters use a recurrent inadmissible increase of the time delay td ¼ nT0/2. InP N1 aðkÞxðn  kÞ and the output signal procedure, the filter equations become simpler yðnÞ ¼ k¼0 depends on the input samples without any dependence on the previous output samples. The class of these digital filters is known as finite-impulse response (FIR) filters, which are always stable and have a response length for an unit impulse d(n) equal to N (Rabinder, Gold [47]). So the standard filtering task consists in the determination of the weight a(k) in the Pm1coefficients aðkÞxk after z-transtransfer function of a digital filter, which has the form HðzÞ ¼ k¼0 formation. At present a large number methods for a(k) determination are known from which averaging windows have been used (Cappellini et al. [48], Gutnikov [49]). One needs here to mention only the problems with the use of such digital filters. The main principal ones are: (1) there appears an inadmissible huge time delay wN/2 because of the large size of the measured sampling signal that is necessary for filtering in the ultralow frequency region and (2) any random impulse in the filtering input noise produces P a the catastrophic filter reaction, which represents the filter’s own pulse response yðnÞ ¼ m1 k¼0 aðkÞdðn  kÞ . This long and strong input disturbance leads to false feedback output control in automated Cz and correspondingly to the crash of the complete digital control system. 5.1. Digital Filtering in Real Time using a Kalman-Bucy Filter The more convenient approach is based on a parametric model of the measured signal and the calculated estimation of its state variable b xðkÞ using the measured sensor signal y(k). Such a filtering algorithm, proposed by Kalman corresponds to the Gauss-Markov recurring estimation procedure. It reduces in general, to a recurring estimation byðkÞ with a minimum root-

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

60

mean-square error (see book Brammer, Ziffling [50]). The output weight sensor signal can be represented by the mathematical model: xðk þ 1Þ ¼ AxðkÞ þ hðkÞ yðkÞ ¼ CxðkÞ þ uðkÞ

ð133Þ

Where k  0, x(k) e is the state vector (inaccessible to direct observation) a noise model, y(k) e is the measured vector, h(k) and u(k) are the vectors of perturbation and measuring error, which are the random processes of white noise type. The state vector in our case can be _ € _ formulated by three component xðkÞ ¼ ½PðkÞ; PðkÞ; PðkÞthe weight value P(k), the ratePðkÞ € and the acceleration PðkÞ of variable weight. Note here, that the scalar value of the measuring error u(k) can be represented as white noise since one can shift its spectral density by two orders to the high-frequency area compared to the frequency cut-off [46]. Here we will only discuss the three-dimensional Kalman filter. This filter allows one to get all the necessary components of the input deviations, which are needed for PID feedback control. In general, such a Kalman filter is represented by the equation:   b xðk þ 1Þ ¼ Ab xðkÞ þ KðkÞ yðkÞ  CAb xðkÞ ð134Þ Where b xðk þ 1Þ e is the current restored value of a state variable, K(k) is the matrix of the gain constant, which can be calculated by equation [50]: 1  KðkÞ ¼ PðkÞCT CPðkÞCT þ R

ð135Þ

Here R is the noise u(k) covariance matrix with zero average main Mu (k) ¼ 0 and P(k) is the covariance matrix of the measuring error which is determined by the equation: PðkÞ ¼ AfPðk  1Þ  Kðk  1ÞCPðk  1ÞgAT

ð136Þ

Before calculating using the recurring eq. (133) it is necessary to consider an initial state. The convergence algorithm is derived for the case of the assumption about the zero average Mx(k) ¼ 0 distribution and the infinity of dispersion value P1 ð0Þ ¼ 0. In practice it is easy to determine the real average value MxðkÞs0 and the noise dispersion P1 ð0Þs0 before Kalman filtering. Also it is possible to refresh at every filtering cycle the knowledge about the noise stochastic characteristics, due to the growing conditions and the current crystal weight. _ x2 and b x3 e are the rate PðkÞ and If b x1 is the weight estimation P(t) and correspondingly b € acceleration PðkÞ of variable weight estimations, so eq. (133) can be represented as: 1 0 b x1 ðk þ 1Þ 1 T @b x2 ðk þ 1Þ A ¼ @ 0 1 b x3 ðk þ 1Þ 0 0 0 1 @0 0 0

1 0 1 10 b K1 ðkÞ x1 ðkÞ 0:5T 2 x2 ðkÞ A þ @ K2 ðkÞ AfyðkÞ  ð 1 0 0 Þg T A@ b b x3 ðkÞ K3 ðkÞ 1 1 10 2 b x1 ðkÞ T 0:5T x2 ðkÞ A 1 T A@ b b x3 ðkÞ 0 1

ð137Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

61

Fig. 37. Data processing during Kalman filtering. (a) - Sensor data for 200 sec measurements with sampling time 1 sec; (b) - 1-recovery of current weight P(t) using a sampling time of 10 sec, 2-the background noise; (c) - Comparison of the _ real filtered weight 1 with that from program 2; (d) - recovery of the weight changes rate PðtÞ.

Where T is the discrete time measurement, K1 ðkÞ ¼ p11 ðp11 þ RÞ1 ; K2 ðkÞ ¼ p21 ðp11 þ RÞ1 ; and K3 ðkÞ ¼ p31 ðp11 þ RÞ1 . The components pij for the measuring error covariance matrix can be calculated using eq. (135) with the initial state Pð0Þ ¼ Að D 000D=20000D=400ÞAT , where D ¼ R e is the current “refreshed ” dispersion sensor signal. Within one filtering cycle it is possible to accept the noise dispersion value as constant. The components p22(0) and p33(0) of the covariance matrix were chosen equal to D/20 and D/400 corresponding to every filtering cycle. These values provide the maximum filtering effectiveness for the three-dimensional Kalman filter. Fig. 37 shows only one filtering cycle at the end of which the current crystal weight and its rate are obtained free from noise; then they are used for immediate feedback control. The curve 2 visually shows the effect of transposition of the frequencies with increasing discrete sampling time [42]. Due to this stroboscopic effect there is a transport of noise to lowfrequency area which no longer allows any more the use of the approach of white noise. In Fig. 37c we compare the program and restored weights at the moment when the velocities _ PðkÞ are the same but the absolute weights P(k) are different, the typical situation in PID feedback control at the initial (or final) growth stage. It is possible to cancel (“forget”) this integral error and recalculate and then select the appropriate program for the current restored cone weights. The above approach allows one to achieve more than 98% noise suppression (Fig. 38). 6. Modelling for Liquid Encapsulated Czochralski growth Developing an appropriate control system for Liquid Encapsulated Czochralski (LEC) growth requires consideration of several additional factors. Since the melt meniscus and

62

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 38. Bar graph of the filtered weight amplitudes in relation to the amplitude of the noisy weight (as a percentage). The following filtering methods were used: 1-Lancosh filter; 2 - Dolf-Chebushev filter with the integer coefficients; 3 - Dolf-Chebushev filter; 4 ethe fast digital inverse Fourier transform (FFT); 5- the flat-apical window; 6-simple averaging; 7-digital filter based on 25 points; 8 - Blackman-Herris filter; 9–Hann window; 10-Hamming window; 11-the two-dimensional Kalman filter; 12-the three- dimensional Kalman filter;

crystallization front are surrounded in the LEC process by liquid encapsulant B2O3 , the melt capillary constant is equal to a2 ¼ 2sLF =Drg where sLF -is the surface tension of the meltencapsulant meniscus surface, Dr ¼ rl  r is the difference in densities of the melt and the encapsulant. The encapsulant also changes the heat exchange around the meniscus and interface area, and an Archimedean force which acts on the crystal and the meniscus to change the measured weight is present. However, similar to the basic Cz process the main effect on the crystallization front and the shape of a growing crystal in the LEC process are melt temperature changes at the base of the meniscus. Several different conditions that occur during LEC growth are analyzed in this section and are shown in Figs. 39 and 40. An additional variant is when the crystal is completely immersed within the liquid encapsulant eFull Encapsulated Czochralski (or FEC).

Fig. 39. Scheme for Liquid Encapsulated Czochralski (LEC) growth using weight control. Changes at the crystallization front contain the derivative of the weight signal.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

63

Fig. 40. In LEC the output weight signal consists in two parts Y ¼ Y1 þ Y2 and this current measurement is retained in the computer memory. The second term Y2(t  td) of the current output signal Y2(t) depends on the crystal diameter cross-section at the upper bound of the encapsulant and includes the time delay td.

Also there are changes to the heat exchange coefficient but similarly to the basic Cz process the main effects on the crystallization front and the shape of a growing crystal in the LEC method are the melt temperature changes at the meniscus base. Linearization of the three conservation equations enables one to formulate the crystallization dynamics in the LEC and FEC methods using a non-autonomous (time-dependent) system which is given by eq. (60) (see Section 3). So the LEC and FEC LOM have the second order and in general the LEC method in closed state can _ be described in the standard matrix form by the equation XðtÞ ¼ AðtÞXðtÞ þ BðtÞðUðtÞ  UðtÞÞ, where UðtÞ- is a perturbation vector, the input vector U(t) ¼ [u1,u2]T also has two components for the control vector u1 ¼ DT(t) and u2 ¼ DV0(t) in the T- and V- (required) feedback channels. However, evaluation of the coefficient matrices A(t) and B(t) is different for the LEC case and is derived in this section. The principal feature of the LEC (not FEC) method consists in weight control in the presence of the encapsulant on the melt surface. 6.1. Weight Control Measurement Analysis for LEC Evaporation of the encapsulant or main material loss are presumed not to occur. In contrast to the basic CZ process the weight signal in LEC includes a time dependent Archimedes force. It is evident that to take into account this component in the current output signal it is necessary to know exactly the “crystal’s prehistory”. This means that it is necessary to know the crystal shape which is immersed under the encapsulant as well as the current meniscus volume. It is important to note, that the melt meniscus does not have a “memory” about its previous states and geometry and is specified by the current crystal diameter and the slope of its lateral surface at the crystallization front. For steady crystallization r0 ¼ const, the lateral surface should be _ were vertical b(t) ¼ 0 and the crystallization rate is equal in simplest case to VC0 ¼ V0  HðtÞ, _ HðtÞ ¼ rs r02 V0 =ðrl R2  rs r02 Þ. Since the production of A3B5 semiconductors (GaAs, InP, GaP) uses rigid crucibles, one does not need to take into consideration the crystallization rate variation due deformation of the crucibles R ¼ const (see Section 3). The output signal of the weight sensor in the case of steady crystallization (see Fig. 39) is determined as: Z t Z t 2 2 _ ð138Þ HðtÞdt  Fa ðtÞ FðtÞ ¼ prs g r ðtÞVC ðtÞdt þ gmðtÞ  Fa ðtÞ ¼ prs gR 0

0

64

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Here m(t) - is the current meniscus mass and: Fa ðtÞ ¼ grf mf ðtÞ=rs þ grf mðtÞ=rl ¼ Fam ðtÞ þ Fam ðtÞ

ð139Þ

is the variable Archimedes force. The crystal section mf (t) and the melt meniscus are immersed in the liquid encapsulant with a density rf. The derivative of the output signal provides the information about the current crystal diameter variations at the crystallization front. Differentiation eq. (137) gives: _ ¼ prs gr 2 ðtÞVC ðtÞ þ gmðtÞ _  F_ a ðtÞ FðtÞ

ð140Þ

For steady cylindrical crystal growth the meniscus mass and its volume are constant m(t) ¼ const and also Fa(t) ¼ const, r(t) ¼ const, VðtÞ ¼ VC0 ¼ const. Therefore the equation _ eq. (139) simplifys to FðtÞ ¼ F_ 0 ¼ prs gr02 VC0 ¼ prl gR2 H_ 0 . Any small changes in the steady crystallization conditions due to melt temperature or pulling rate variations causes the local crystal shape variation drðtÞs0 and accordingly the change of the meniscus height dhðtÞs0. So the output weighing signal to second order accuracy is equal to:   _ ¼ prs g r2 þ 2r0 drðtÞ ðVC0 þ dVC ðtÞÞ þ gdmðtÞ _  dF_ a ðtÞ ð141Þ F_ 0 þ dFðtÞ 0 The output equation describes the deviations from the known program value F_ 0 and can be written as:   _ _ ¼ g dmðtÞ _ þ dmðtÞ  dF_ a ðtÞ ð142Þ dFðtÞ _  dF_ a ðtÞ ¼ gdMðtÞ _ since the rate value of the crystal mass variation is equal to dmðtÞ ¼ prs ðr02 dVC ðtÞ þ 2r0 drðtÞVC0 Þ with second order accuracy. The derivative of the Archimedes _ _ f ðtÞ=rs Þ þ grf ðdmðtÞ=r _ force variation dF_ a ðtÞ ¼ grf ðdMðtÞ=r l Þ ¼ grf ðdm l Þ may be represented by the pair components dF_ a ðtÞ ¼ dF_ am ðtÞ þ dF_ am ðtÞ. Now we use an important simplification in the calculations: assume, that rf =rs yrf =rl ¼ h sincers yrl . Finally we obtain the output equation for the LEC method using weight control as the very simple output equation which is similar to that for the basic Cz method: _ ¼ gð1  hÞdMðtÞ _ dFðtÞ

ð143Þ

This eq. (142) describes the output signal of a weighing sensor for the current variations of the crystal diameter and the meniscus volume at the crystallization front, immersed in an encapsulant (see Fig. 39). The Archimedes force is time dependent due to the changes in the crystal (dmf /rs) and meniscus (dm/rl) volumes within the liquid encapsulant. Now let us take into account the Archimedes force variation dFam for any current time t ¼ t0, when the perturbed crystal leaves the encapsulant. In this case the rate of crystal mass variation (immersed in the encapsulant) can by written as: 0 1 Z t0 Z t0 td _ 0 Þdt  _ 0 ÞdtA dmðt dmðt ð144Þ dFam ðt0 Þ ¼ ghdmf ¼ gh@ 0

0

Here the value td determines the moment in time before the current time t0 > td, when the crystal section is generated, which coincides with current time t0 at the upper encapsulant boundary (Fig. 39). The first integral in eq. (143) describes the complete crystal mass variation at the time t0 and the second integral describes the variation of the crystal mass immersed in the liquid encapsulant. The variation of meniscus volume at the current time t0 determines the second part

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

65

dFam ðt0 Þ ¼ ghdmðt0 Þ of the full Archimedes force variation dFa and does not depend on the crystal’s prehistory. After differentiation of the described values dFam(t0) and dFam(t0) and substitution of the results in eq. (4.141) we obtain finally the measurement equation at time (t0) as: _ 0 Þ ¼ gð1  hÞdMðt _ 0 Þ þ ghdmðt _ 0  td Þ dFðt

ð145Þ

or using the mass balance in time t ¼ t0  td , i.e.: _ 0  td Þ ¼ dmðt _ 0  td Þ þ dmðt dMðt _ 0  td Þ This eq. (144) can be represented otherwise as: _ 0 Þ ¼ gð1  hÞdMðt _ 0 Þ  ghdMðt _ 0  td Þ  ghdmðt dFðt _ 0  td Þ

ð146Þ

In standard matrix form the output equation for the LEC method can be written as: YðtÞ ¼ Y1 ðtÞ þ Y2 ðt  td Þ ¼ C1 XðtÞ þ D1 UðtÞ þ C2 Xðt  td Þ þ D2 Uðt  td Þ ð147Þ Here t ¼ t0 is the current time of measurement. The matrices C1(t) ¼ [c11, c12] and D1(t) ¼ [d11, d12] look like the time dependent matrices for the Cz model: C1 ¼ ðc111 c112 Þ ¼ PK1 ½ðAHr _ þ AH_ r_ a11 þ AH_ h_ a21 ÞjðAH_ r_ a12 þ AH_ h_ a22 Þ D1 ¼ ðd111 d112 Þ ¼ PK1 ½ðAH_ r_ b11 þ AH_ h_ b21 ÞjðAH_ r_ b12 þ AH_ h_ b22 þ AHV _ Þ

ð148Þ

but here PK1 ¼ prl R2 ð1  hÞ and the capillary constant in LEC is equal to a2 ¼ 2sLF =Drg. The delayed component in the output eq. (146) is determined by the matrices C2(t) ¼ [c21, c22] and D2(t) ¼ [d21, d22] with the coefficients: c211 ¼ pghrs r0 ð2V0  r0 Arh Þ c212 ¼ pghrs r02 Ahh d211 ¼ pghrs r02 AhT d212 ¼ 0

ð149Þ

So, in contrast to the basic Cz method, there appears in the LEC method a delayed component to the weighing signal Y2(t  td)which depends on the height of the “additional” liquid-gas boundary at the upper encapsulant boundary) Fig. 40. This component “encrypts” the information about the current crystal and the meniscus variations at the crystallization front. Note here, that the measurement (output) eq. (146) is correct not only for cylindrical crystal growth, but also for the process of pulling of an initial- or a final crystal cone (the monitoring mode for feedback). The output equation in this case remains the same, but its matrices become time dependent. Taking into account distortion of the crystallization front and the variations of its volume U(r(t)) it is possible for them to be similarly represented as in Section 2.2. The derived eq. (60) given in Section 3 and (146) represent the inertia-free linear mathematical model (LOM) for Liquid Encapsulated Czochralski. Below we take into account the dynamics of the melt temperature response u1 ¼ DT as the result of changes of heating power DN. 6.2. Trajectory and the time delay calculation Fig. 41 shows the scheme for pulling an initial cone in the LEC method. The algorithm which is used as a program task (trajectory for the monitoring mode of a controllable feedback

66

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 41. (a)-Scheme for pulling an initial cone in the LEC process. (b)- The surface of trajectories Y(b(r(t))), as a function of the predetermined initial cone shape b(r). One of them (bold curve) will represent the program trajectory for a closed-loop system.

system) is based on the use of a recurrence procedure. Prior to the crystal immerging from the encapsulant the following integral equations are correct. At any time the weighing sensor determines the weight: PðtÞ ¼ g½mðtÞ þ mðtÞ  Fa ðtÞ þ const The Archimedes force determined as:   dFa ðtÞ ¼ gh mf ðtÞ þ mðtÞ

ð150Þ

ð151Þ

where mf(t) - is the mass of crystal, which is covered with encapsulant. From these equations and the mass balance MðtÞ ¼ mf ðtÞ þ mðtÞ it follows: PðtÞ ¼ gð1  hÞMðtÞ

ð152Þ

If the current total crystal mass is more then the mass covered with encapsulant m(t) > mf(t), i.e. the crystal emerges from the encapsulant it follows:    ð153Þ PðtÞ ¼ g mðtÞ þ mðtÞ  h mf ðtÞ þ mðtÞ ¼ gð1  hÞMðtÞ þ ghmðt  td Þ since MðtÞ ¼ mðtÞ þ mðtÞ and mf ðtÞ ¼ mðtÞ  mðt  td Þ: The eqs. (151) and (152) describe the real weight signal structure P(t). The time derivative determines the current state of the dynamical system for weighing control. For the program task value it is necessary to rewrite eqs. (151) and (152) as b_ ¼ gð1  hÞMðtÞ b_ PðtÞ

ð154Þ

b_ ¼ gð1  hÞMðtÞ b_ b_  td ðtÞ PðtÞ þ ghm½t

ð155Þ

Here and below we denote the program values using the symbol “cap” L. Now we will calculate the program task values for obtaining the crystal with a predetermined shape using the appropriate function b(r), which describes the slope angle to the crystal profile. The following b_ b_ _ _ ¼ prl R2 HðtÞ ¼ V0  where VC(t) and HðtÞ equations mðtÞ ¼ prs r 2 ðtÞVC ðtÞ and MðtÞ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

67

VC  ðh0r þ h0b b0r ÞVc tgb are derived in Section 4.2. Also it is necessary to store the results in the computer memory after every calculation step. The total number of stored trajectory (Fig. 41b) values n depends on the choice of the discrimination time T0 for digital feedback. One accepts as a first step i ¼ 1 and k ¼ 1 (definition of counter k is given below): r1 ¼ rin ¼ const at the next step rm, m > i using the function b(r) define bi ¼ b(ri), then one can calculate the current meniscus height hi(ri, bi) using eq. (4) (Section 2.1) and the current crystallization rate VCi ¼ eq. (108) (Section 4.2). So the current crystal radius on the i-th step VC ðri ; bi Þ usingP P is equal to riþ1 ¼ rin þ T0 Vci tgðbi Þ and the crystal current length is equal to li ¼ T0 VCi . Correi i b_ ¼ b_ , spondingly the time derivative of the current crystallization mass is equal to M prl R2 H i i Pb 2 b ¼T _ the total mass M and the current meniscus mass is m b ¼ pr r h M i

i

0

i

i

l i

i

bi  m bi ¼ M ðri ; bi Þ þ prl a2 cosð3 þ bi Þ. Finally the current crystal mass m bi is calculated and is stored in the computer memory. Determination of the time delay td(t), which depends on the variable encapsulant height Hf (t)is fundamentally important for the LEC method. Defining this time td(t) and assuming the constant encapsulant mass in the crucible Mf ¼ const. To determine Hf (t) one needs to look at the total volume under the R t encapsulant in the crucible with a radius R ¼ const. The full crystal length is equal to lðtÞ ¼ 0 VC ðtÞdt) at the time t. Before the final emergence of the seed crystal from under encapsulant the correct relation is:   ð156Þ pR2 Hf ðtÞ þ hðtÞ ¼ Wm ðtÞ þ Wmf ðtÞ þ Wf 0 þ Ws where Wf 0 ¼ Mf =rf ¼ const is the encapsulant volume, Wm ¼ m bðtÞ=rl is the meniscus volume, b f =rs is the volume of the crystal covered by encapsulant and Ws ¼ prs2 ½Hf  lðtÞ is Wmf ¼ m the volume of the seed crystal rs, covered by encapsulant. Taking into account eqs. (149) and (150) we obtain the mass balance solution:    b 1 m bðtÞ rl R2 hðtÞ MðtÞ 2   1  þ W Hf ðtÞ ¼  2 þ  pr lðtÞ  2 ð157Þ f0 s 2 rs rl rs p R  rs R  rs2 At the moment in time, when the top of the initial cone emerges into the gas environment from under the encapsulant l(t) ¼ Hf (t) and Ws ¼ 0, so:    b 1 MðtÞ m bðtÞ rl Hf ðtÞ ¼ 1  þ Hf 0 þ hðtÞ þ ð158Þ pR2 rs rl rs where Hf 0 ¼ Mf =prf R2 . And when the growing crystal appears from under the encapsulant it is b f ðt  td Þ: necessary to use the stored values m    b  mf ðt  td Þ m 1 MðtÞ bðtÞ rl Hf ðtÞ ¼ 1 þ Hf 0 þ hðtÞ þ ð159Þ rs pR2 rs rl It is possible to represent the derived formulas (156)-(158) as HF ðtÞ ¼ Hf ðtÞ þ hðtÞ ¼ Hf 0 ðtÞ þ f ðM; m; m; tÞ. From physical considerations it is clear, that the encapsulant thickness can increase or decrease depending on the cone shape for materials having a density relation rl =rs > 1, for example for GaAs. By way of illustration we show in Fig. 42 the calculated encapsulant thickness for two different initial crystal cone shapes.

68

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 42. The calculated encapsulant thickness Hf (r(l )) can increase or decrease depending on the initial cone shapes. The initial encapsulant thickness Hf 0 ¼ 1.2 cm.

In order to determine the time delay value td (t) it is necessary to take into consideration the time of Rseeding and to calculate and remember the crystal length l(t) and the crystal mass t b _ b ðtÞ ¼ 0 mðtÞdt. If the next calculation step performs the inequality l(t) > Hf (t) this means m that the crystal is beginning to emerge from under the encapsulant. It is necessary to restore td (t) using stored values and eq. (128). The correct determination of the value td (t) consists in controlling the equality lðtÞ  lðt  td Þ ¼ Hf ðtÞ. Now to continue the algorithm of the program task for the LEC method: the calculated thickness of the encapsulant Hfi (t) is determined by eq. (156) for Hfi > li, or by the eq. (157) in the case that Hfi ¼ li and finally using eq. (5.158) when Hfi < li. So the program trajectory will be calculated from the equation: b_ ¼ gð1  hÞM b_ P i i

ð160Þ

when the initial cone is completely covered by encapsulant and by the equation: b_ ¼ gð1  hÞM b_ þ ghm b_ ik ; where k ¼ td =T0 P i i

ð161Þ

when the top of the initial crystal cone emerges from under the encapsulant. It is evident, these eqs. (159) and (160) are the discrete analogue of the eqs. (153) and (154) since a discrete feedback process is used to calculate the above values in a discrete time t ¼ iT0. The counter k in eq. (160) determines the stored values which are needed to determine the variable time delay td. 6.3. PID control for LEC processes This section examines whether a single channel PID controller works well for the LEC process. When the initial cone is completely covered by encapsulant, feedback control design in

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

69

the LEC process is similar to the basic Cz design (see Section 4). In this case, the weight output eq. (144) for LEC does not have the delayed component, and the matrix notation looks similar to the output equation in basic Cz (see Section 2.4) YðtÞ ¼ Y1 ðtÞ ¼ C1 ðtÞXðtÞ þ D1 ðtÞUðtÞ

ð162Þ

They differ only by the predetermined coefficient (1  h), where h ¼ rf =rs xrf =rl . The value U (t)in eq. (161) is the vector of the input perturbations in the T- and V-cannels. Due to the similarity of the LEC and the Cz transfer functions in this case one can refer for the detail to Section 2.5. As well it in the LEC model is necessary to take into account an inertia of the melt temperature response to controlling heating power changes DN (t). Similarly it is necessary to add to the equations any additional transfer functions such as first or second order melt temperature responses for example. The difficulties in using a simple PID control after the crystal cone emerges from the encapsulant comes from the fact that a weight sensor contains simultaneously the information about the current as well as the previous variations of the crystal diameter and meniscus volume: _ ¼ gð1  hÞMðtÞ _ _  td Þ ð163Þ PðtÞ þ ghmðt Due to this second delayed component there is a possible control error, which is absent in the basic CZ method. If time t is the current moment in time, so that at the previous time t  td, when the growing crystal was located under the encapsulant (m ¼ mf) the following values were measured and stored in the computer memory: Pðt  td Þ ¼ gð1  hÞMðt  td Þ ¼ gh½mðt  td Þ þ mðt  td Þ   _  td Þ ¼ gð1  hÞMðt _  td Þ ¼ gh mðt _  td Þ þ mðt Pðt _  td Þ

ð164Þ ð165Þ

dPðt  td Þ ¼ gð1  hÞdMðt  td Þ ¼ gh½dmðt  td Þ þ dmðt  td Þ   _  td Þ ¼ gð1  hÞdMðt _  td Þ ¼ gh dmðt _  td Þ þ dmðt dPðt _  td Þ

ð166Þ ð167Þ

The last eqs. (165) and (166) describe the output variations which can be used for feedback control at the time t  td. These variations e are the deviations of the current values from the predetermined programmed values: b f ðt  td Þ ¼ dmf ðt  td Þ mf ðt  td Þ  m

ð168Þ

b_ f ðt  td Þ ¼ dm b_ f ðt  td Þ m_ f ðt  td Þ  m

ð169Þ

Taking into account eq. (168) the current measured signal of a weight sensor (162)) is:   b_  td Þ þ dm_ f ðt  td Þ _ ¼ gð1  hÞMðt _  td Þ þ gh mðt ð170Þ PðtÞ From this equation it follows that even if the time delay td (t) is known exactly, it is impossible determine the current output control object variations. Let’s assume, at the current time t no perturbations take place near to the crystallization front and any deviations from the program _ ¼ 0, i.e. they are completely absent. As usual, the desired program value value dMðtÞ ¼ dMðtÞ b_ b_ b_  td ðtÞÞ can be subtracted from the current measured value PðtÞ _ PðtÞ ¼ gð1  hÞMðtÞ þ mðt b_ _ _  PðtÞ (eq. (169)) which can result in a nonzero value dPðtÞ ¼ PðtÞ ¼ dm_ ðt  t Þs0. So, in the f

d

70

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

LEC process it is possible to have incorrect feedback control using a simple PID controller which is not caused by external disturbances close to the crystallization front. Unfortunately it is impossible to determine the value dm_ f ðt  td Þ in the stored output signal, described by eq. (169). The inverse of the trajectory calculation, i.e. restoring the crystal and meniscus shapes using only the output weight signal, does not provide an accurate solution. The difficulty of this task consists in the necessity to register all the input control and perturbations which produce measured deviations. While determining the desired trajectory, it is sufficient to use the predetermined function b (r) and the mass balance (see Section 4.2), for recovering the crystal and meniscus shapes, it is necessary to take in to account the heat balance at the crystallization front, i.e. to use the complete description of the crystallization dynamics. That suggests, however, that in general it is possible to solve this task by using an observer (see Section 6.4). In practice the main errors in the LEC method, especially at the time of growth of the initial cone, are determined by an incorrect time delay td (t) determination. Indeed, before the crystal cone emerges from under the encapsulant, the value of the input control is determined by small deviations (165) and (166). After that, when the crystal cone begins to emerge from under the encapsulant incorrect feedback control depends first of all, on the accuracy of the determination of the time delay td (t). Formally this inaccuracy can be estimated as: b_  td Þ þ dm_ f ðt  td Þ _ þ DPðtÞ _ ¼ gð1  hÞdMðt _  td Þ þ ghDmðt dPðtÞ

ð171Þ

b_  td Þ  mðt b_  td  dtd Þ is the error of the program value due to incorrect b_  td Þ ¼ mðt where mðt b_  td Þ can choice of the time delay td ðtÞ  dtd. It is evident, that the error in the program value Dmðt b_ f ðt  td Þ  dm b_ f ðt  td  dtd Þ. be much more than the current measured error dm_ f ðt  td Þxdm Physically this means, that the cross-section of a real cone at the upper encapsulant boundary differs from the calculated cone cross-section. So, the main feature of the LEC method consists in the superposition of a measured output signal with information about the current controlled object state _ dMðtÞ _ 0 Þ þ hdmðt _ 0 Þ=g ¼ ð1  hÞdMðt _ 0  td Þ YðtÞhdFðt

ð172Þ

_  td Þ state at the simultaneously with the additional information about the subsystem dmðt previous time. It is important to note, that the time delay in general depends on the measurement time td ðti Þhtdi . For the PID controller to work correctly it is necessary in the output signal YðtÞ ¼ Y1 ðtÞ þ Y2 ðt  tdi Þ to separate the component Y1 (t) depending on the current crystal state XðtÞ ¼ ½drðtÞ; dhðtÞT . Initially consider the case for pulling a cylindrical crystal, as is shown in Fig. 43. Consider as a first approximation that the value Hf (t) varies insignificantly, that is - the time delaytdi xHf =VC0 ¼ const, were i- is the number of measurement moments. The simplest approach consists in the determination of a current mass deviation directly from measured output signal Y (t) (see eq. (170)): 1 h _ _  td Þ ð173Þ YðtÞ þ dmðt dMðtÞ ¼ 1h 1h _  td Þ represents the error in the determination In this equation, the unmeasured value dmðt of the current deviation rate of crystal and meniscus mass from the desired value. However, one can store the measurements and use them to improve the approximation accuracy of the current output state at the crystallization front.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

71

Fig. 43. During cylindrical crystallization, to a first approximation, one uses the supposition that the time delay is constant tdi xHf =VC0 ¼ const.

It is possible to represent the information for the current (t ¼ ti) measurement by the following equation: _ _ td Þhdmðt _ _ td Þ ¼ ð1hÞdMðtÞhd Mðt _ td Þ ð174Þ YðtÞ ¼ ð1hÞdMðtÞþhd mðt The stored measurement at the previous step ti1 ¼ ti  td (see Fig. 43) e is accordingly to: _  td Þ þ hdmðt _  2td Þ Yðt  td Þ ¼ ð1  hÞdMðt _  td Þ  hdMðt _  2td Þ  hdmðt ¼ ð1  hÞdMðt _  2td Þ

ð175Þ

And at the next previous step ti2 ¼ ti  2td - correspondingly to: _  2td Þ þ hdmðt _  3td Þ Yðt  2td Þ ¼ ð1  hÞdMðt _  2td Þ  hdMðt _  3td Þ  hdmðt ¼ ð1  hÞdMðt _  3td Þ

ð176Þ

and so on. So the following combination can be made from the current Y (t) and the stored Y (t  td) output signals of the weight sensor (Fig. 44): YðtÞ 

h h2 _ _ dMðt  2td Þ  hdmðt _  td Þ þ Yðt  td Þ ¼ ð1  hÞdMðtÞ 1h 1h h2 dmðt _  2td Þ þ 1h

ð177Þ

_  2td Þ and dmðt _  2td Þ for materials with h ¼ rf =rs xrf =rL < It is possible to neglect dMðt 0:4 (for GaAs hx0:3). The next number step in this combination does not change the final result: YðtÞ 

h h2 _ Yðt  2td Þ ¼ ð1  hÞdMðtÞ Yðt  td Þ þ  hdmðt _  td Þ þ Oð3Þ 1h 1h

ð178Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

72

_ Fig. 44. An accurate approximation of YðtÞhdF=g(eq. (170)) based on eq. (176). In approximation 1 all terms for the current Y (t) and the previous value of the weight signals (stored) Y (t  td) Y (t  2td) were used. A step-wise dT ¼ þ4C for GaAs (see Table 3) was used.

So, only two sensor signals- the current Y (t) and the stored Y (t  td) allows one to estimate _ the deviations at the crystallization front dMðtÞ as:  1 h h _ ð179Þ YðtÞ  Yðt  td Þ  dmðt _  td Þ dMðtÞ ¼ 1h 1h 1h In this combination the value dmðt _  td Þ represents the determination of the error of the current deviation at the front from the program value. It is possible to demonstrate, that eq. (177) _ approximates to the value dMðtÞ better than eq. (171). Substituting in this equation the value _ 0  td Þ ¼ Y (t  td) from eq. (172) and taking into account the mass balance dMðt _ 0  td Þ þ dmðt _ 0  td Þ we obtain: dmðt _ dMðtÞ ¼

1 h h2 _  td Þ  _ YðtÞ þ dmðt 2 dmðt  td Þ 1h 1h ð1  hÞ

ð180Þ

_  td Þ if we only use two signals, the current The approximation error is reduced in 0:24dmðt Y(t) and the stored Y(t  td) (Fig. 45). In a similar way it is possible to describe weight control during the growth of an initial cone. Fig. 46 shows the initial stage of crystal growth in an LEC process. Before the crystal cone begins to emerge from under the encapsulant the output signal _ _ YðtÞ ¼ ð1  hÞdMðtÞ ¼ ð1  hÞðdmðtÞ þ dmðtÞÞ _ allows one to save in memory the value: _ dmðtÞ ¼

1 YðtÞ  dmðtÞ _ 1h

ð181Þ

which will be used later at time t ¼ t þ td . At time t ¼ t0 þ td ðtÞ after pulling the crystal into the _ _  td ðtÞÞ contains the stored gas phase the measured weight signal YðtÞ ¼ ð1  hÞdMðtÞ þ hdmðt

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

73

Fig. 45. In approximation 2 two output signals were used e the current Y(t) and one stored value Y(t  td). The other parameters are the same as in Fig. 44. h _ signal (179). As result for the initial LEC stage we get YðtÞ  1h Yðt  td ðtÞÞ ¼ ð1  hÞdMðtÞ  hdmðt _  td ðtÞÞ which is completely identical to (176), but here the time delay td (t) is in principle the variable value.

6.3.1. Transfer processes with closed loop PID control of the LEC process If feedback control in LEC is based on using current deviations from the program at the _ so it is possible to use known solutions for autocrystallization front Y1 ðtÞ ¼ ð1  hÞdMðtÞ, mated Cz. When the LEC crystal leaves the encapsulant it is possible to use the values Y (t) and Y _  td ðtÞÞ, so for the real approximation(t  td). But since nothing is known about the value hdmðt h _ Yðt  td ðtÞÞ. This sum was 2 the accuracy is defined by the equation ð1  hÞdMðtÞ ¼ YðtÞ  1h calculated and is shown as Approximation 3 in Fig. 47. So it is necessary to compare the quality of PID control based on Approximation-3 using only the current measured output signalY(t)¼Y1(t)þY2(ttd). Consider the scalar step

Fig. 46. Structure of output signal using weight control in LEC during the growth of the initial cone.

74

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 47. Approximation 3 for the current deviation at the crystallization front as the sum ð1  hÞ _ dMðtÞ ¼ YðtÞ  h=1  hYðt  td ðtÞÞ of the current Y(t) and the stored Y(ttd(t)) weight measurements. The parameters and crystal material are the same as in Fig. 44.

perturbation of the melt temperature dT ¼ þ 4C for the growth of a GaAs crystal with the parameters listed in Table 4. The calculated transfer processes for open (without feedback control) state LEC are shown in Fig. 48. One can evaluate the performance of closed loop control for such a system using a PD controller with a single T-control channel. The weight deviations for the program values (I-component) are not used in our calculations. Here, the current “measured” output signal is Y Table 4 Parameters for GaAs using the LEC crystallization process Technological parameters (well-known) Seed crystal radius Crystal radius Radius of crucible Angle of crystal cone Crystal pulling rate

(m) (m) (m) (grad) (m/sec)

Rs R0 Rcr B V0

4.5*0.001 40*0.001 100*0.001 50 4*(0.001/3600)

Material parameters Thermal conductivity solid phase Thermal conductivity liquid phase Latent heat Growth angle Capillary constant Density of liquid phase Density of solid phase Density of encapsulate

(W/mK) (W/mK) (J/kg) (grad) (m) (kg/m3) (kg/m3) (kg/m3)

ls ll L E0 A rl rs rf

17.12 7.2 728700 15 4.2/1000 5.71*0.001/(0.01)^3 5.17*0.001/(0.01)^3 1.75*0.001/(0.01)^3

Technological parameters (rough estimate for GaAs crystallization process) Melt overheating (K) DT Temperature gradient in crystal (K/m) Gs

4 1500

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

75

Fig. 48. Calculations for a LEC GaAs crystal (r ¼ 0.95r0) with transfer processes for open state after a scalar step-wise perturbation dT ¼ þ4C. (a) e As a result the crystal irreversibly transfers from A into a new stationary state B; (b) e components of the weight output signalY(t)¼Y1(t)þY2(ttd)with a time delay td ¼ 22min. (e) e the crystallization rate variations cannot be more than dVc100% (see Section 3.3.3)

(t)¼Y1(t)þY2(ttd). So, the digital feedback control consists in the calculation of the control _ dTc ¼ kp YðtÞ þ kd YðtÞ at every sampling timeT0, which can compensate the consequence of the step-wise melt temperature increase in dT ¼ þ 4C. Also used is a PD digital regulator, which uses two weight signals instead of a single current weight output signal. The difference will consist in the use of approximation-3 (see Fig. 47) based on the current Y(t) and the saved Y (ttd) signals (and its derivatives). The initial parameters and controller coefficients kp and kd remain the same. Fig. 49 shows the transfer processes calculated for such a system. The proposed analysis demonstrates the possibility of successful working control with a simple single channel parametric regulator which only uses the T-channel with a PD and a simple structure (without integrator) for crystal diameter stabilization with the LEC method. Moreover it is shown that it is possible to use only a signal Y(t)¼Y1(t)þY2(ttd) for current deviations from a program. There is some improvement in control quality, when the approxi_ mation ð1  hÞdMðtÞ ¼ YðtÞ  ðh=ð1  hÞÞYðt  td ðtÞÞ is used instead of the single signal Y (t) (see Fig. 44) which can not compensate for the difficulties obtained in control. So, in general, diameter control in the LEC method does not differ from the use of a PID regulator

76

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

_ Fig. 49. Closed loop response, calculated for an LEC process with PD-single channel T-control dTc ¼ kp YðtÞ þ kd YðtÞ _ for a GaAs crystal. (a) e use of the approximation ð1  hÞdMðtÞ ¼ YðtÞ  h=1  hYðt  td ðtÞÞ for current deviations at the crystallization front; (b) e using only the current measured output signal Y(t)¼Y1(t)þY2(ttd) as the P- and its derivative (the second derivative of weight) as the D-component of the PD-controller. The discrete controlling time T0 ¼ 60sec and the parameters and perturbation are as in Fig. 38.

with the Cz method. Poor choice of PID coefficients in the feedback controller can result in crystal diameter oscillations as shown in Fig. 50. These oscillations appear for large values of the P- and D-control parameters when the remaining technological parameters for the GaAs growth process remained the same.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

77

Fig. 50. (a)- Response oscillations for closed loop LEC with single T-channel PD feedback control with enlarged coefficients. The calculations used: time delay td ¼ 22min, material GaAs, crystal radius r ¼ 0.95r0 (see also Table 4). (b) e photo of LEC GaAs crystal with the incorrect feedback set up for a single T- channel PID controller.

The problem for the correct set up for a PID controller will be considerably complicated due to an unsuccessful construction of the heating zone and a correspondingly bad temperature distribution in the melt and in the thermal zone for crystal growth . Fig. 51 a shows a similar growth zone which had a temperature distortion near to the middle of the melt level in the crucible. For any small melt overcooling/overheating dTw0 C under the crystallization front the LEC close loop response dr(t) becomes strongly oscillatory (Fig. 51b) in response to the temperature u1s0 (and also the pulling rate u2s0) input. Therefore the PID controller, which was set up

Fig. 51. Example of a bad construction for a growth zone in LEC. (a)-In addition to the temperature T(L) distortion (here near to the bottom of the crucible) also: an unsuitable crucible bottom shape and irregular crucible heating due to the step-by step deterioration of the graphite support properties; (b)- calculated envelope for a GaP crystal response dr(t) (as a function of the melt overheating/overcooling dTs0) for the input u1 ¼ dT ¼ þ5C. Here the parameters: rcr ¼ 030min, V0 ¼ 20mm/h, Gs ¼ 150C/cm (see text).

78

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 52. A visual representation for a digital control system which allows one to calculate and show the response r(t)¼ r0þdr(t) for the LEC process. The calculated crystal shape in an open- (right in image (a)) and closed loop state (right in image (b)) corresponds to the temperature distribution, shown in Fig. 51a. Also the required (cylindrical) crystal is shown left as an image. (c) e A GaP crystal, obtained under similar conditions with normal feedback control during pulling of its first half length and bad control for the remaining lower length. (d) The work operator for a high pressure puller using the digital crystal diameter system.

correctly and worked successfully for the upper part of crucible loses its stability when the solidified front falls below the temperature distortion (Fig. 51a). In practice this effect appears suddenly as a large and rapid variation in crystal diameter (Fig. 52c, e) and for the operator (Fig. 52d), without any “visible” reason . Very often the operator undertakes an unhappy attempt to suppress such an instability and often with a resulting separation of the crystal from the melt. This example illustrates the extraordinary importance of the solution for solving the problem defining the variable parameters of dynamical systems and the necessity for adaptive feedback control. This last problem is currently not solved for practical basic Cz and LEC crystal growth.

6.4. Observers and state control in LEC processes The discrete LOM model for the LEC process looks like:  Xðk þ 1Þ ¼ FXðkÞ þ fððUðkÞ  UðkÞÞÞ YðkÞ ¼ CXðkÞ þ DðUðkÞ  UðkÞÞ þ C1 Xðk  nÞ þ D1 ðUðk  nÞ  Uðk  nÞÞ

ð182Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

79

Fig. 53. The block-scheme for LOM in the state space for the LEC process using weight control. It is similar to the Cz block-scheme (see Section 4.4.1 and Fig. 31) before the appearance of a LEC crystal from under the encapsulant. And when the LEC crystal leaves the upper bound of the encapsulant the output signal changes and consist of two parts: the current Y1(k) and with delay on n tracts Y2(kn) (so Y(t)¼Y1(t)þY2(ttd)).

Where T0 ¼ const, n ¼ td(t)/T0 the variable discrete time delay, FðT0 Þ ¼ expðAT0 Þ ¼ PN n 1 n¼0 ðAT0 Þ =n!, f(T0)¼(FI) AB, I-is the unit matrix 2x2. (discussed in Section 2.6.) . The diagram in state space for the LOM (179) for the LEC process in a closed loop state is represented in Fig. 53. Our controlled object is completely observable and controllable in the case of double channel T- and V-control. It allows one to use the results, obtained for basic Cz in Section 4.4.1 and build the equivalent Luenberger observer: 8 b þ1Þ ¼ FXðkÞþfUðkÞþHDeðkÞ b < Xðk b mÞþD1 Uðk mÞ b b ð183Þ YðkÞ ¼ CXðkÞþDUðkÞþC1 Xðk : b DeðkÞ ¼ YðkÞ YðkÞ b b Here XðkÞ and YðkÞ e are the calculated values of the state variables and an output, H e is the matrix which defines the upper bound of the encapsulants observer convergence rate. It is important, that the time delay for the observer m is not equal to the real time delay n but is a variable task parameter. The error for state reconstruction will be determined as:   b  mÞ þ HD1 ðUðk  nÞ b þ 1Þ ¼ ðF  HCÞDXðkÞ þ HC1 Xðk  nÞ  Xðk DXðk  Uðk  mÞÞ

ð184Þ

This error depends both from the common rate of convergence lim ðF  HCÞDXðkÞ ¼ 0 k>N and the error in the selection of an appropriate value m compared to the real discrete time delay n. The observer convergence rate is determined by the roots of the characteristic equation: detðzI  F þ HCÞ ¼ 0

ð185Þ

Corresponding to the separation theorem the closed state dynamics of the controlled object is determined by the set of the characteristic equation det(zIFþfK)det(zIFþHC)¼0 poles

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

80

b [44]. So, the required observer matrix H and the state controller UðkÞ ¼ KXðkÞ matrix K can be determined separately. The procedure for calculating this matrix is described in detail in Section 4.4.1. The estimations of the state variables can be calculated if one knows the input U(k)and output Y(k) signals of a controlled object as: b  mÞ  HD1 ððUðk  mÞ b þ 1Þ ¼ ðF  HCÞDXðkÞ b þ ðf  HDÞUðkÞ  HC1 Xðk DXðk þ HYðkÞÞ

ð186Þ

So the estimation error consist in two parts De(k)¼De1(k)þDe2(km). Minimization of error b  mÞ. De2(km) is based on correct choice in the memory of the stored values U(km) and Xðk The quality of the work of the observer depends on an informal operation, the selection of the time delay m. 6.4.1. Determination of the time delay in LEC td(t) The total encapsulant height (above the melt level in a crucible) is determined as: HF ðtÞ ¼ Hf ðtÞ þ hðtÞ

ð187Þ

for a crystal growing with radius r(t) and meniscus height h(t), value Hf(t) e is the thickness of the encapsulant layer measured from the three-phase position (Figs. 40 and 41). In the case of the growth of an initial cone Hf(t) was defined in Section 6.2 (eqs. (156) and (157)). For stationary crystallization with 2r0 ¼ const it is possible to estimate the value td0 ¼ H/Vc0 to a first approximation, Hf0 e is the encapsulant thickness. This stationary encapsulant height Hf0 is determined from the solution of the equation set: HF0 ¼ Hf 0 þ h0 2 2 r1 f Mf ¼ pR HF0  pr0 h0 þ pa r0 cosð3Þ  pr0 HF0

ð188Þ

as the value:   Hf 0 ¼ Mf =prf R2  r0  a2 r0 cosð3Þ=R2  r0  h0

ð189Þ

More accurately the time delay td(t)¼td0þdtd(t) after taking into consideration the variations dHf(t) due to variations of the meniscus volume dWm(t), crystal volume dWC(t) and also the crystallization rate dWC(t). If at any time t ¼ 0 perturbations occur then after the time interval t ¼ t0the point at the “perturbed ” crystal surface is shifted relatively to the crystallization front a distance: Z t0 dVC ðtÞdt ð190Þ lðtÞ ¼ VC0 t0 þ 0

and when this cross-section reaches the upper encapsulant level the value Hfbecomes equal to: Z t0 Hf ðtÞ ¼ Hf 0 þ dHf ðtÞ ¼ VC0 td þ dVC ðtÞdt ð191Þ 0

Here td ¼ t0 e is the time when this crystal cross-section reaches the upper encapsulant level.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

81

In general, when the current time interval exceeds the time delay value t0>td(t0), then the eq. (188) can be rewritten as: Z t0 dVC ðtÞdt ð192Þ Hf ðt0 Þ ¼ Hf 0 þ dHf ðtÞ ¼ VC0 td ðt0 Þ þ t0 td

So the time delay td(t)¼td0þdtd(t) is determined by the current encapsulant thickness Hf(t0), which depends on the volume variations of the meniscus and the crystal, also from the crystallization rate variations. From the eq. (189) it follows: 2 3 Z t0 1 dVC ðtÞdt5VC0 ¼ td ðt0 Þ þ dtd ðt0 Þ ð193Þ td ðt0 Þ ¼ 4Hf 0 þ dHf ðtÞ  t0 td

From the current values dHf(t0) and dtd(t0) it is possible to determine from an analysis of the volume balance, which is situated between the melt and the encapsulant levels (Fig. 41): W ¼ Wm þ Wmf þ Wf 0

ð194Þ

where Wm ¼ m/rl is the current meniscus volume, Wmf ¼ mf /rS e is the volume of crystal under the upper encapsulant level and Wf0 ¼ Mf /rf e is the constant encapsulant volume. Since   W ¼ W0 þ dWðt0 Þ ¼ W0 þ pR2 dHf ðt0 Þ þ dhðt0 Þ  mr drðt0 Þ þ mh dhðt0 Þ Wm ¼ Wm0 þ dWm ðt0 Þ ¼ Wm0 þ r1 l Z t0 Z t0 ð195Þ _ _ Wmf ¼ r1 dmðtÞdt mðtÞdt ¼ pr02 VC0 td þ r1 S S t0 td

t0 td

so taking into account the relation W0¼Wf0þWmf0þpr20VC0td0, from eqs. (190)e(192) we obtain: 2 3 Z t0 drðtÞdt  pR2 dhðt0 Þ5 ð196Þ dHf ðt0 Þ ¼ l1 4dWm ðt0 Þ þ 2pr0 VC0 t0 td

2 dtd ðt0 Þ ¼ l2 4dWm ðt0 Þ2pR2

Z

t0 t0 td

dVC ðtÞdt 2pR2 dhðt0 Þþr1 S

3

Z

t0 t0 td

_ 0 Þdt5 dmðt

ð197Þ

where l1¼1/p(R2r0), l2 ¼ l2/VC0 andRthe partial derivatives m0r and m0h were derived in Section _ and t0 dhðtÞ _ ¼ dhðt _ 0 Þ  dhðt _ 0  td Þ so this eq. (194) is 2.1. SincedVc ¼ dV0  dH_  dh, t0 td possible to be formulated as: 22 2 dtd ðt0 Þ¼l2 44r1 l dmðt0 Þ2pR

þr1 l

Z

Z

t0

t0 td

dV0 ðtÞdt2pR2 dhðt0 Þþr1 S

Z

t0

t0 td

_ 0 Þdt dmðt

3 t0

t0 td

5 _ dMðtÞdt

ð198Þ

82

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Rt Here the value dV0 is the control in the V-channel. Since t00 td dmðtÞdt _ ¼ dmðt0 Þ  dmðt0  td Þ _ ¼ 0, so after simplification rSxrl in eq. (195) we finally obtain: and dm_ þ dm_ þ dM 2 3 Z t0 2 dV0 ðtÞdt  2pR2 dhðt0  td Þ5 ð199Þ dtd ðt0 Þ ¼ l2 4r1 l dmðt0  td Þ  2pR t0 td

Using the state variables x1 ¼ dr(t), x2 ¼ dh(t) and control in the V-channel u2 ¼ dV0in (196) we obtain: Z t0 2 u2 ðtÞdt ð200Þ dtd ðt0 Þ ¼ C1 x1 ðt0  td Þ þ C2 x2 ðt0  td Þ  pR t0 td

Where C1 ¼ l2 mr =rl and C1 ¼ l2 ðm0h =rl  pR2 Þ. So, the variation of the time delay is determined by the value of the state variable X(t0td) and the output control in the V-channel u2 ¼ dV0, if it is used in a feedback control and determined by the stored state variables. Let a discrete control system use the sampling timeT0. Suppose, at the i-th moment in time t ¼ iT0 it is known that m ¼ i  k are the values of the state variables, where k ¼ 1,.,(i1). It is necessary determine the value td(t)¼td0þdtd(t). That means it is necessary to “remember” (store) and then calculate the value: k X u2 ð201Þ dtd ðiÞ ¼ C1b x1 ðmÞ þ C2b x2 ðmÞ  pR2 T0 i

and simultaneously calculate the divisible of the estimated time: nðkÞ ¼ ðtd þ dtd ðiÞÞ=T0

ð202Þ

In general searching for k continues while fulfilling the condition: n  m < ðn þ 1Þ

ð203Þ

It is obvious, that it is necessary to carry out this searching near to the values k ¼ i  m0, where m0 ¼ t0/T0. More precisely the state variables can then by determined by the equation < y; r >¼ 0 Where x e is the real rest in eq. (199) and k ¼ n. If in a digital control system it does not use the V e control channel, the value u2 ¼ 0 in eq. (198). 7. Weight based control The development of high quality systems for controlling crystal diameter and for the identification of controlled object parameters requires additional analysis of the adequacy of sensors and their real performance. However, in many papers, control of the CZ and LEC process is based on the assumption of ideal sensitivity of the sensors without taking into account their real characteristics. This section analyzes this problem for the weight control method, as well as the problem of sensor interaction with the controlled object and complete measurement of the dynamical process. 7.1. Weight control Weight control is the most commonly used indeed universal practice for crystal growth since it has minimal effect and any disturbance is symmetrical with the heating zone. Weight sensor

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

83

design and any interaction with the dynamical system is therefore very important. Among measuring devices the simple tensometric sensor can be referred to as a first order object. This implies a static error during a variable load measurement. Since the output value of such a sensor depends on the absolute magnitude of the deformation of its string element, so the variable load involves some lag from the final string position. Such a final position is used to calibrate a tensometric sensor in a static situation. More interesting are the second order measuring devices. Their advantage consists in the capability of adjusting the dynamics of a precision device’s characteristics. This section analyzes the interaction of a self-compensated balance (SCB) sensor with a growing crystal (Cz) or a crucible containing a melt and is based on articles and patents [51e54]. The general requirements of the weighing devices are: (1) a guarantee of the independence of its precision with a load of variable value, (2) the measurement results must be independent of the position of the loading point and (3) the mechanical part of the self-compensated balance must be in indifferent equilibrium in the open state. The principal features of the SCB sensor mechanics are shown in Fig. 54. The SCB sensor uses a short elastic hinged (“spring-like”) strip and a mechanical adder with three or more levers so that there is only vertical movement from the remaining one degree of freedom. The short elastic strip allows the hinge to guarantee the achievement of requirements (1), (2) and additionally achieves maximum sensitivity and measurement linearity. This is possible, since the restoring momentum M after deformation of the short elastic strip hinge (Fig. 54a) depends weakly on the loading value (Vasil’ev at al. [55]):   1=2 ð204Þ M ¼ ðPEJÞ cthðxÞ ¼ EJl1 1 þ x2=3 þ x4=45 þ . where P is the loading value, E e the modulus of elasticity of the strip of hinge material, l e the qffiffiffiffiffiffiffiffiffiffiffiffi e is length of short elastic strip hinge, J e moment of inertia of its cross-section, x ¼ l Pcosð4Þ EJ the reduced strip hinge length and 4 e the hinge deviation angle from an initial position. For x2, the dependence (201) becomes practically linear for a large weighing diapason. The third condition e the indifference to the equilibrium of the self-compensated part of the mechanical balance is possible due to the compensation of the elastic strip hinge’s inflexibility

Fig. 54. (a) Scheme for the self compensated balance (SCB weight sensor). (b) Photograph of the sensor’s mechanism that fulfills the scheme’s concept.

84

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 55. Scheme which explains the working and mechanics of computation of the SCB weight sensor; (a) the restoring moment of the elastic rotating joint M is compensated by shifting of loads on the heights H1 and H2 (b)-variant of a device for rotating from rod 5 with a scheme of a mechanical hammer and an elastic rotating joint; (c) e alternative of this device, which uses friction and centering between the element 6 and 7. This device eliminates breakage of the central rod in the case that the rotating crystal freezes in a crucible.

due to the shifting of the lever’s centre of gravity Z in position H1 relative to the hinge jamming point N(see Fig. 55a). In the SCB device shown in Fig. 54b this moment of force of the compensation is created by the load arrangement on the levers. The variable strip hinge’s inflexibility is dependent on the variable load. It is also compensated after the additional shifting of the lever’s centre of gravity in position H2. This position depends on the position of the instantaneous rotation axis Q which acts on the resultant forceP0þF. The detailed calculation for an optimum elastic strip hinge involves the shape and heights H1 and H2 and is based on the relation M(P)¼PpH1þ(P0þF )H2, where Pp e is the levers weight including the compensating loads and P0 e is the constant load of the receptacle. Figs. 56 and 57 shows the model of the “upper” SCB weight sensors.

Fig. 56. (a)- 3d scheme of an “upper” SCB weight sensor; (b)- photo of the sensors operating mechanism.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

85

Fig. 57. The block-scheme of an SCB sensor in a closed loop state

The SCB sensor includes a set of interrelating elements e a magneto electric force element, an electronic feedback control system and a lever position sensor (Fig. 57). Every unit can be represented by a transfer function. So the transfer function for the photoelectric lever position sensor W2 at the zero initial balanced state is equal to W2(s)¼ DI2(s)/Dx(s)¼K2/(1þT2s), where Dx-is the movement or shifting lever,DI-the electrical photoelectric sensor output, T2 e its time constant and K2 is the transfer ratio. The magneto electric force device may be represented by the transfer function W3(s)¼DF(s)/DI3(s)¼K3, were DI3 is the input electric current in the solenoid, DF- the force developed in the magneto electric force device and K3 ¼ gBln is its transfer ratio, depending on the magnetic induction B in the air-gap, ln the full winding length of solenoid and g e the lever reduction coefficient. In order to determine the main transfer functionW1(s)¼DI3(s)/DP(s), which connects an input (steady or variable load) with an output (the electric current or voltage in solenoids), it is necessary to analyze the device dynamics when it deviates from the initial balanced equilibrium state. 7.2. Constant load Let us consider measuring a constant load of mass ms by a SCB sensor moving vertically upward with a constant rate V0 (Fig. 54a). Initially, when the sensor is situated in the equi€ ¼ 0 < so for the accuracy of the residual restoring moment U it librium state 4 ¼ 4_ ¼ 4 follows that: m0 g  gF ¼ 0

ð205Þ

Here m0 is the constant mass of the sensor elements controlling the measurements. It is possible to derive the set of motion equations for the sensor using d’Alembert’s equilibrium principle: 8 Pn > > < 1 F=n ¼ 0 ðms þ m0 Þg  dtd ½ðms þ m0 Þ  Nz ¼ 0 ð206Þ > d > _  U4=n ¼ 0 pðN : z =nÞ  pgðF=nÞcos4  dtððI0 =nÞ4Þ Here g- the ratio of a sensor lever length (the reduction coefficient), n- the lever numbers, Vs-the moving rate of a constant load, Nz is the supporting force, I0 is the constant moment of the levers inertia and 4 e the small levers deviation angle. The first equation in (202) is the force equilibrium, the second and last equation is the consequence of the equality of the moments of force relative to n lever rotating axis On and its projection in the coordinate system

86

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

(z,r). The solution of these equations and the constraint equation Vs ¼ p4_ gives the equation of motion:   € HU4 ð207Þ mgp ¼ I þ ms p2 4 where m ¼ ms þ m0 and I ¼ I0 þ m0p2 is the constant which characterizes the sensor’s inertial properties. The presence of a residual moment U is the main cause of sensor faults. A positive residual moment þU accounts for oscillation in the dynamic system. The result obtained determines the lack of perfection of the strain sensors. From eq. (203) it is clear that information about value of constant load P ¼ mg, gives the movement acceleration of the levers and the first derivative (rate of movement), defines the variable load since the lever shifting x w 4. The growing crystal (or crucible) weight measurement is more complicated due to the interaction of a pair of dynamic systems since the load is in principle variable. 7.3. Variable load; crucible weight The system’s initial equilibrium state will be considered for the case of pulling a cylindrical crystal of constant diameter 2r0, so the crystallization rate is equal Vc0 ¼ V0  H_ 0 and the derivative of the mass of the growing crystal m_ ¼ prs r02 Vc0 ¼ const (Fig. 58a). As a result of € ¼ 0 of sensor, connected with a crucible, mass balance 4i , the equilibrium state 4 ¼ 4_ ¼ 4 having the initial melt mass M¼prlR2H0, is described as: _ 0 VM0 ¼ 0 ðM0 þ m0 Þg  gF  M

ð208Þ

_ 0 =2prl R2 is the stationary movement rate of the crucible with its melt where VM0 ¼ 0:5H_ 0 ¼ M and H_ 0 is the initial melt level (see Fig. 58a).

Fig. 58. Weighing scheme. (a) Measurement of the crucible mass M(t); SCB sensor is static (fixed); the crystal moves with a constant rate V0. (b) Measurement of the crystal mass m(t) and meniscus mass m(t); the SCB sensor and the crystal move with a constant rate V0.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

87

Assuming that in time dt the load changes like: DM ¼ dM þ dM

ð209Þ

_ 0 dt- the stationary decrease of liquid in the crucible and dM¼prlR2dH is the here dM ¼ M _ variation of melt mass associated with the weighing process 4s0 and the input disturbances of the crystallization process U(t)¼[dT(t),dV0(t)]Ts0. A set of equilibrium equations according to the D’ Alembert principle for variable mass ms(t) are: 8 Pn F=n ¼ 0 1 > >  < d d ðms þ m0 Þg  ðms Vs Þ  m0 p4_  Nz ¼ 0 ð210Þ dt dt > > d  : pNz  pgFcos4  I0 4_  U4 ¼ 0 dt Here ms(t)¼M0þdM(t)þdM(t), Vs(t)¼VM0þdVM(t) e the movement rate of the center of _ þ p4_ is the movement rate variation for the crucible mass of the load and dVM ðtÞ ¼ 0:5dHðtÞ plus the centre mass of the melt. The equation of motion, derived from eq. (206) can be rewritten finally as:     _ 0 4_  U4 ð211Þ €  0:5pM _ 0 þ prl R2 VM0 dH_ ¼ I þ 0:5p2 M 4 € þ p2 M pgDM  0:5pMdH with second order accuracy. Here I ¼ I0þpm0 and M(t)¼M0þdM(t)þdM(t), and it is assumed _ M ¼ Oðd2 Þ. dMdV 7.4. Variable load; crystal weighing The weighing scheme is shown in Fig. 58b. The instantaneous equilibrium state for the € ¼ 0, coupled to the crystal and moving vertically at a constant rate V0, is sensor is4 ¼ 4_ ¼ 4 described by the following equation: _ m0 ¼ 0 ms ð0Þg  gF þ Fm0  mV

ð212Þ

where ms0 ¼ m0 þ m(0) is the total initial mass of the load including the crystal mass m(0) and the constant mass of the crystal holder, taking part in the accelerated movement at weighing. The value Fm0¼m0g¼pgr20h0þa2grlcos(30) is the instantaneous equilibrium meniscus weight. At the state of stationary crystallization m_ 0 ðtÞ ¼ 0 and Fm0 ¼ 0. Here h0 is the stationary meniscus height, a¼(2sLG/rg)1/2 is the Laplace constant, m_ ¼ prs r02 V0 ¼ const and Vm0 ¼ 0.5Vc0 is the stationary displacement rate of the crystal’s centre of mass (see below). Assume that in time dt the load changes by: ms ðtÞ ¼ ms ð0Þ þ dmðtÞ þ dmðtÞ

ð213Þ

Fm ðtÞ ¼ Fm0 þ dFm

ð214Þ

_ is the stationary change of the crystal mass, dm(t) and dFm =g ¼ m0r drðtÞ þ Here dmðtÞ ¼ mdt e are the mass variations of the crystal and meniscus, associated with the weighing _ process 4s0 and input control or perturbations in the T- and V- channels of the controlled object. Now it is necessary to define the displacement rate of the crystal’s centre of mass, having an arbitrary shape (non cylindrical lateral surface). Again, assume that the crystallization front is

m0h dhðtÞ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

88

flat or does not changes its curvature and, taking into account a sensor response 4(t) as a result of the crystal mass increase Dms(t)¼dms(t)þdms(t). Now it is necessary to take into account the relations (Fig. 58b): _ Z_ 1 ðtÞ ¼ p4ðtÞ _ _  dhðtÞ Z_ 2 ðtÞ ¼ V0  H_ 0 þ dVðtÞ þ dHðtÞ _ Z_ 3 ðtÞ ¼ V0  H_ 0 þ dVðtÞ  dHðtÞ Z_ 4 ðtÞ ¼ V0 þ dVðtÞ _  p4ðtÞ _lðtÞ ¼ Vc0 þ dVc ðtÞ ¼ V0  H_ 0 þ dVðtÞ  d HðtÞ _  dhðtÞ _

ð215Þ

The displacement rate of the crystal’s center of mass Z_ m ðtÞ, assume that the crystal has a cylindrical shape on average and drdrmax¼const for ti˛[0,t]. At the measurement time the crystal radius is r(t)¼rþdr(t). In general the position of the centre of mass is determined as: Z Z2 o n  2 Z Z2 2 pr r ðtÞ  z ðtÞ þ 2 zdrðz; tÞdz z 0:5r 0 0 s 2 1 zr 2 ðz; tÞdz Z1 Z1 ð216Þ ¼ Zm ðtÞ ¼ Z Z2 ms ðtÞ 2 r ðz; tÞdz Z1

Here, dm(t)<
z1

where Z*˛[Z1,Z2] belongs to the interval l(t)¼Z2(t)Z1(t) equal to the crystal length. From eqs. (210), (212) and (213) it follows that:   _ þ p4ðtÞ _ ð1  drmax =r0 Þ ð219Þ Zm ðtÞ ¼ 0:5 V0  H_ 0 þ dVðtÞ  dhðtÞ And finally the displacement rate of the crystal’s centre of mass is obtained with second order accuracy is:   _  dHðtÞ _ þ p4ðtÞ _ Vm ðtÞ ¼ Z_ m þ dZ_ m ðtÞ ¼ 0:5VC0 þ 0:5 dVðtÞ  dhðtÞ ð220Þ where Z_ m ¼ Vm0 ¼ 0:5VC0 is the stationary part and dZ_ m ðtÞ ¼ dVm ðtÞ ¼ 0:5½dVðtÞ _  dHðtÞ _ _ dhðtÞ þ p4ðtÞ is the variation of a load centre of mass displacement. A set of equilibrium equations according to the D’ Alembert principle for a variable growing crystal mass ms(t) are: 8 Pn > F=n ¼ 0 > < 1  d d  ðms þ m0 Þg þ Fm ðms Vs Þ  m0 V0 þ dV  p4_  Nz ¼ 0 > dt dt > : pN  pgFcos4  d ðI 4Þ _  U4 ¼ 0 z

dt

0

ð221Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

89

Substituting here the equations (208), (209) and (215), the equation of motion is obtained with an accuracy of second order dm_ s dVm ¼ Oðd2 Þ in the form: _ m ðtÞ  0:5pVC0 dmðtÞ _  pm0 dV_ m ðtÞ ¼ I€ 4ðtÞ  U4ðtÞ pgDms ðtÞ  pms ðtÞdV_ m ðtÞ  pmdV ð222Þ or finally as:

  € þ 0:5pmd _ þ 0:5pms ðtÞd€hðtÞ þ 0:5p m_ þ VC0 m0 d€hðtÞ _ HðtÞ gpDms ðtÞ þ 0:5pms ðtÞdHðtÞ h _  0:5pmdVðtÞ _ þ o:5pVC0 m0r d_r ðtÞ  p½0:5ms ðtÞ þ m0 dVðtÞ   2 2 € ðtÞ þ 0:5p m_ 4ðtÞ _  U4ðtÞ ¼ I þ 0:5p ms ðtÞ 4 ð223Þ

where Dms ðtÞ ¼ dmðtÞ þ dmðtÞ þ dmðtÞ and VðtÞ ¼ V0 þ dVðtÞ. Again the dynamical system is of second order with variable parameters. Note there is a difference in crucible and crystal the weighing dynamics. More complicated transfer dynamics for crystal weighing are caused by the moment of force as a result accelerated crystal and meniscus centre of mass movement. In the case of crucible weighing, this moment of force is compensated for by a supporting force (the upper puller rod), which is hard connected with the crystal. Note here, that direct crystal weighing (the "upper" measurement scheme) excludes the errors due to the crucible levitation with induction heating. 7.5. Feedback control sensor design The differential eqs. (207) and (218) allow one to define the transfer function of the SCB sensor in the open- and closed states. The coefficients are variable and change during weight measurements. But it is possible to consider the assumption of “frozen coefficients” assuming ms ðt þ dtÞ ¼ ms ðtÞ if the time required for measurement is small dt/N. For further analysis of the weighing dynamics, it is convenient to use the Laplace transform of eqs. (207) and (218), assuming the initial conditions to be zero. The Laplace transformation of eq. (207) gives: 4ðsÞ ¼ WM ðsÞDMðsÞ þ WHM ðsÞdHðsÞ

ð224Þ

And in the case of weighing of a growing crystal from eq. (218) one obtains: 4ðsÞ ¼ WM ðsÞDms ðsÞ þ WHm ðsÞdHðsÞ þ Whm ðsÞdhðsÞ þ Wrm ðsÞdrðsÞ þ WVm ðsÞdVðsÞ The transfer functions derived from eqs. (219) and (220) looks like:

ð225Þ

_ 0  Ug ¼ gp=DM ðsÞ WM ðsÞ ¼ gp=fs2 ðI þ Mð0Þp2 Þ þ sp2 M    _ WHM ðsÞ ¼  0:5s2 pMð0Þ þ s pMð0Þ þ prl R2 VM0 DM ðsÞ 2 2 2 Wm ðsÞ ¼ gp=fs ðI þ 0:5ms ð0Þp Þ þ 0:5sp m_  Ug ¼ gp=Dm ðsÞ Wrm ðsÞ ¼ 0:5spVC0 m0r =Dm ðsÞ _ WHm ðsÞ ¼ f0:5s pms ð0Þ þ 0:5spmg=D m ðsÞ Whm ðsÞ ¼ f0:5s2 pms ð0Þ þ 0:5spðm_  VC0 mh 0Þg=Dm ðsÞ _ WVm ðsÞ ¼ fsð0:5pms ð0Þ þ m0 Þ þ 0:5pmg=D m ðsÞ 2

ð226Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

90

Since the LOM of the Cz growing system (see Section 3) taking into account the weighing sensor may be written as: 8 < d_r ¼ Arr dr þ Arh dh _ dh_ ¼ Ahr dr þ Ahh dh  dH_ þ AhT dT þ ðdV0  p4Þ ð227Þ : _ _ dH ¼ AHr dr þ AHr_ d_r þ AHh_ dh_ þ AHV ðdV0  p4Þ So, according to the Laplace transformation of eq. (222) one can obtain the transfer functions of the controlled object taking into account the presence of the weight sensor. The solution of the equations set: 3 2 3 3 2 Arh 0 drðsÞ s  Arr 0 4 s  Ahh s 5 4 dhðsÞ 5 ¼ 4 AhT dTðsÞ þ dVðsÞ  sp4ðsÞ 5 Ahr sAHr_  AHr sAHh_ s AHV ðdVðsÞ  sp4ðsÞÞ dHðsÞ 2

ð228Þ

yields: drðsÞ ¼ WrT dT þ WrV ½dVðsÞ  sp4ðsÞ dhðsÞ ¼ WhT dT þ WhV ½dVðsÞ  sp4ðsÞ dHðsÞ ¼ WHT dT þ WHV ½dVðsÞ  sp4ðsÞ

ð229Þ

The transfer functions in eq. (224) were obtained above (see Sections 1.4 and 1.6). It is obvious that the joint resolution eqs. (221) and (224) describe the process dynamics for the growing crystal or crucible weighing including the presence of noise and an external input U(t). It is convenient to determine the measurement error of the SCB sensor in the closed loop state as: eðtÞ ¼ gDms ðtÞ  FðtÞ

ð230Þ

Determination of the sensor feedback circuit structure involves finding the polynomial order of the transfer functions WfB ðsÞ ¼ KðsÞ=ZðsÞ and coupling the angular deviation of the sensor levels to the compensating force FðsÞ=g ¼ WfB ðsÞ4ðsÞ. Without loss of generality, one can assume transfer functions for the position sensor of the levers and force element for the SCB to be equal to unity. Analyzing the closed-loop dynamics for the most commonly used practice when the object input is the temperature only, dTs0 and dV¼0, weighing the crucible with the melt is determined by: 8 < 4ðsÞ ¼ WM ðsÞ½dMðsÞ þ dMðsÞ þ WHM ðsÞdHðsÞ dHðsÞ ¼ WHT ðsÞdTðsÞ  spWHV ðsÞ4ðsÞ ð231Þ : eðsÞ ¼ DMðsÞ  WfB ðsÞ4ðsÞ One obtains:

    eðsÞ ¼ e1 ðsÞ þ e2 ðsÞ ¼ ð1 þ spWHM WHV ÞDMðsÞ  WfBM WHM WHT dTðsÞ VM

ð232Þ

where VM ¼ 1 þ spWHM WHV þ WfBM WM . Precise sensor operation requires a zero steady state error lim eðtÞ ¼ 0, caused by the melt temperature variation dTs0 and correspondingly the t/N sensor input DMðtÞ ¼ dMðtÞ þ dMðtÞ. Since for stationary crystallization dM(t) changes linearly it is necessary that the minimum compensation by the sensor feedback for standard power input tuðtÞ ¼ 1=s2 , where u(s) is the unit step.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

91

The simplest controller based on a SCB sensor is a PID compensator WfBM ðsÞ ¼ ðf0 s2 þ f1 s þ f2 Þ=s. The simple verification, based on eq. (232) and the final-value theorem shows, that the static error lim eðtÞ ¼ lim e1 ðsÞ þ lim e2 ðsÞ ¼ 0 in case of precise compensation of t/N s/0 s/0 a residual moment U¼0 and a linear disturbances dM(t). Since the transfer processes in a selfcompensated balance are fast (see below), the linearity of dM(t) is a good approximation. The sensor response for a more complicated input perturbation, for example, dTðtÞ ¼ T0 tn uðtÞ=n!, n¼1,2,. or the presence of a residual moment Us0 has some static errors which are easily estimated as: lim e1 ðsÞ ¼ ðArr Ahh þ Arh AHr  Arh Ahr ÞU=gpf2 _ 0 Arh AhT AHr =g lim e2 ðsÞ ¼ M

s/0 s/0

In general, improving the measurement accuracy is possible by using the PII2 for SCB feedback. But in practice, using an additional integral component in the feedback leads to instability. Optimization of a SCB sensor feedback R t is based on the calculation of its transfer process and using the integral quality criterion J ¼ 0 e2 ðtÞdt ¼ min. The technique for such calculations was described in [34]. In Fig. 59 are shown the results of a typical transfer process, based on eq. (227) and the parameters of an SCB sensor shown in Fig. 56a,b. For visualization of the SCB sensor response we used non-optimum feedback controller parameters and the “slowed” time. The best results correspond to full compensation for a string element residual moment U¼0. In the calculations, enlarged values were used, which in practice are typically smaller than U¼1.105kgm2. Note the remarkable feature of a self compensated balance with a PID controller in feedback. The presence of a small static error (see above) during variable load measurement enables one to

Fig. 59. (a) The block-scheme for crucible weighing (the “lower” measurement scheme), based on eq. (224). (b)- The calculated output response of an SCB sensor after the step-wise melt temperature perturbation dT¼5 C. Initial state is: a stationary (steady) pulling process at a constant rate V0¼10mm/h for a LiNbO3 crystal 80mm in diameter and a mass of 1 kg. The SCB uses a PID controller in feedback. (A)- the balancing error; (B) the output of the displacement sensor of the SCB which allows direct measurement of the rate of loading.

92

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

directly determine the changes of the rate of load. The output electrical signal of the levers displacement sensor may be increased and be used for control of the weighing rate. The larger the weighing rate the greater is the turning angle of the levers and accordingly the output signal. 8. Determination of parameters important for controlled system design High temperature, microwave radiation, vacuum or high pressure in a puller chamber give very restricted access to the growth zone. As a consequence the severe restriction of the experimental possibilities make determination of the controlled object parameters a very difficult task. Here we briefly describe the procedure for determinating the value of the growth angle for opaque and transparent materials, the determination of the capillary constant and the density for high temperature melts. Also we will estimate the dissymmetry of the melt temperature temporal lag in response to heating power changes, the factor of melt temperature expansion and the approach to controlling melt composition. Constant attention to the evolution of the experimental procedures is required. 8.1. Crystallization of melt drops and the growth angle (opaque materials) After separation of a growing crystal from the melt takes place the shrinkage of the remaining liquid film into a small melt drop occurs because the melt wetting of the crystalline surface is incomplete. In these experiments on the crystallization of melt drops it is possible to determine the value of the growth angles at the crystallization front and test some conclusions of the theory. Fig. 60 shows the process of Ge drop crystallization on an isotropic, slightly convex front. In V.V Voronkov’s work [56] he first proposed the use of such cones for the determination of the growth angle. In general, there exist three different cases for the initial location of such drops (and cones) at the crystal front: 1. the cone base (and the melt drop base) is completely located on an isotropic solid surface. 2. the melt drop is completely located on a singular planar facet at the crystallization front (Fig. 62a). 3. the drop base is located both on an isotropic and on a planar singular facet .

Fig. 60. Photos of the step-by-step crystallization of a Ge drop on an isotropic growth front immediately after the crystal is pulled from the melt. The image is inverted with respect to direction of the crystal pulling. The drop formation is completed when the interfacial forces in the plane of the crystallization front reach equilibrium. The vertical component of the summation of the interfacial tensions is compensated by solid phase deformation. (a) eCrystallization front having a curvature with an angle a.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

93

Fig. 61. (a)-Schematic - the cone on an isotropic crystal front (see Fig. 60). (b)-graphical solution of the equilibrium equation gS ðnS Þ ¼ gf ðnf Þ þ gm ðnm Þ for the point L; (c)-position of surfaces in the vicinity of the three-phase line L during drop crystallization, which corresponds to the point A on the Gs figure; (d)-the angular Voronkov diagrams. The changes of the slope angles to the liquid phase at the TPL during melt drop crystallization corresponds to “shifting” from point 1 to 2 on the Gf figure. Accordingly the changes of slope angles to the solid phase vary from point A up to point B on the Gc figure.

In the first case (1) a wonderfully symmetric cone appears after the melt drop crystallizes on an isotropic growth front. This fact allows one to carry out the geometrical construction shown in Fig. 61 [17]. The Wulff Gi figures, where i¼crystal, front and melt, it is possible to represent the situation only qualitatively since the values of the interphase energies and the shape of the polar diagrams as of now are unknown. However, it is possible to take into account some general postulates including the fact that, in this case, singular orientations n0 do not take part in cone shaping. In Fig. 61b is shown the graphical solution of the equilibrium equation gS ðnS Þ ¼ gf ðnf Þ þ gm ðnm Þ for the point L. Using the approach based on a Hoffman-CahnVoronkov construction [5,11] (see also Section ). The centre of the circle Gm “moved” clockwise from the point 1 to 2 along the Wulff figure during crystallization of a melt drop on the apex of the solid cone. The closed triangle (AO1) composed from vectors gi represent the graphical solution of eq. (6). The top of this triangle also “turns” clockwise along Gc from the position A to B. In Fig. 61c is shown the equilibrium positions of interfaces in the vicinity to the three-phase line L, which correspond to the point A in Fig. 61b. The crystallization front near to L should be curved according to the Gibbs-Thomson effect, so there exists the angle 4f to the growth front in the vicinity of L (Fig. 61c). In Fig. 61d are shown the angular diagrams 4c(4m) and 4f(4m) for this case (1). The growth angle is 3¼4m4c and according to Fig. 61d

Fig. 62. (a)- photo of an InSb cone on the singular (111) planar front, (b) e the scheme with its notations used in the calculations

94

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

the growth angle is almost constant. This conclusion follows because in this case (1) the sectors 1e2 and AB on the Wulff figures Gf and Gs can by approximated by circles and Gf should be similar to Gs and moreover the distance between Gf and Gc is practically constant. 8.2. Determination of the growth angle with the assumption of its constancy It is possible to calculate the solid cone shape and compare it with its real shape if one assumes a constant value of the growth angle and the macroscopic front curvature during melt drop crystallization. In such calculations one uses a pair of conservation laws: (i) the mass balance dM ¼ dm þ dm ¼ 0, where dm- is the crystal mass and dm- is the melt mass in the drop and (ii) the growth angle constancy 3 ¼ 4  j ¼ const (see Fig. 62). From the first requirement (i) it follows that rS dVS ¼ rL dV. Since dVS ¼ pr 2 ðzÞdz, finally we obtain: dV=dz ¼ hpr2 ðzÞ, where h ¼ rS =rL . The second requirement (ii) provides the second differential equation dr=dz ¼ ctgðjÞ. For this set of simultaneous equations the initial conditions are the following: the radius of the base of the cone is r(0)¼r0 and its initial volume V(0)¼pr30f(f0)(3þf2(f0)/6, which defines the initial angle f0 (which is identical with the wetting angle in this case). The function f(f)¼(1-cos(f))/sin(f). It is convenient to represent the resulting equations in dimensionless form using the following variables r ¼ r=r0 , z ¼ z=r0 and c¼3V(t)/pr3. After differentiating the last relation dV/dz¼r2c(dr/dz)þr3(dc/dz)/3¼-hr2 and taking into account that dr/dz¼-ctg(j) we obtain the following pair of equations: dc 3 ¼ ðcFðcÞ  hÞ dz r dr ¼ FðcÞ dz

ð233Þ ð234Þ

Here F(c)¼ctg{2arctg(f)-3o} where the initial conditions are: rð0Þ ¼ 1, c(0)¼0.5f(f0) (3þf2(f0). The angle f¼f (c) can be defined as the solution of the cubic equation f3(f)þ3f (f)e2c¼0, which represents the current volume V¼ph(h2þ3r2)/6 by using function f(f). This cubic equation has the unique real solution f(f)¼tg(0.5f)¼[cþ(1þc)1/2]1/3þ[c-(1þc)1/2]1/3 according to the Cardan formula. The solutions of the equations set (228) and (229) depend on the value of two parameters - the growth angle 3 and the initial (wetting) angle fm0 for the crystal material. Fig. 63 shows photos of the solid cones (on the left) and their calculated profile curves (on the right). The real and calculated cone shapes coincide well. This fact confirms the supposition about the growth angle and its constancy in these experiments. It is interesting, that for the determination of the unknown parameters 3 and fm0 it is possible to only use the features of the cone geometry. One such cone feature is the unique

Fig. 63. The cones: (a)-Ge on an isotropic growth front, (b) - Si on an isotropic growth front.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

95

Fig. 64. (a) - Calculated slope angles g(r)¼4(r)-3 of the lateral cone surface of GaAs. The angle at the cone apex g (r/0)¼g0 is constant and does not depend on the initial angle 4o and is determined only by the parameter h¼rS/r and the growth angle 3; (b) e photo of a GaAs cone; (c)- graphical solution of eq. (230) (for a GaAs cone 4*¼69 ).

dependence of the angle at the cone’s apex and on the density ratio h and the growth angle (i.e. material properties of the crystal alone). This angle does not depend on the initial conditions of crystallization. This fact was first mentioned in articles by Antonov [57] and Sanz [58]. Fig. 64 illustrates well such a cone feature. The calculations done above show that the slope angles g (r)¼f(r)-3 of the cone’s lateral surface at the cone apex g ðr/0Þ ¼ g0 are determined only by the choice of the parameters h and 3 and do not depend on the initial angle f0 [17,59]. It is necessary to define an asymptotic solution of eqs. (228) and (229), taking into account, that always 30, f3 and r/0 when z/z0. After differentiation of the drop volume V¼pr3f(f) (3þ f2(f)/6 and after substituting in the derived expression the value dr/dz¼-ctg(j) we obtain:   dV ¼ 0:5DðfÞrz0  0:5rf ðfÞ 1 þ f 2 ðfÞ sin1 ðfÞ40z ¼ h dz

ð235Þ

where D(4)¼f(4)(3þf2(f)). Since for r/0 always 4(r) >0, and so the second term in (230) is zero and can be neglected for the final stage of the crystallization of the melt drop on the top of the cone. Substituting the derivative r0z [-ctg(j) in the remaining part of the expression (230), we obtain the relation 0.5D(f)ctg(j) ¼  h which relates the tangent of the slope angle to the meniscus at the TPL and the slope angle of the crystal’s lateral surface. Taking into account the relation j¼ 4  3 and g¼p/2 j we finally derive the eq. (230) which can be rewritten also as 0.5D(4)ctg(f  3) ¼ h or: 0:5Dð4ÞtgðgÞ ¼ h

ð236Þ

The eq. (231) allows one to determine the growth angle 3 for materials with known h using the measured angle g0 at the cone top: initially from a graphical solution of eq. (231) we determine the value 4* (Fig. 65), which satisfies to the relation D(4*)¼2h/tg(g0) then the growth angle is determined from the expression 3¼4*-(90 -g0). So the measurements give for GaAs the value of the angle at the cone top 2g0¼74 and the graphical solution for this angle gives the value 4*¼69 so the growth angle is 3 ¼16 for GaAs. 8.3. Determination of the growth angle using capillary shaping features The derivation used in the above approach is appropriate and correct for opaque materials, for example semiconductor materials such as Ge, Si etc. For transparent crystals similar to

96

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 65. Graphical solution of eq. (231). For a GaAs cone 4*¼69 . More detailed experiments and singularities for the crystallization of melt drops is described in [17].

Lithium Niobate (Li2NbO3) or sapphire (Al2O3) it is impossible to use separated drops at the crystallization front. In these cases the heat loss is much greater and crystallization of the drop starts on its top. We do not have the chance to obtain a perfect cone at the crystallization front. But it is possible to propose a procedure for the determination of the growth angle 3 for such transparent crystals based on the capillarity properties associated with the crystallization of thin rods by the EFG method (Fig. 66). If for the radius of a growing crystal r/a<1 then the melt capillary constant a ¼ ð2sLG =rl gÞ1=2 , so the meniscus in EFG crystallization will have a large curvature. It is permissible to neglect the hydrostatic pressure for such menisci rl gzðrÞ compared with the Laplace pressure, where z(r) is the meniscus height, measured from its base. For such menisci the Laplace capillarity equation is (in nonedimensional form) rz00 þ z0 ð1 þ z02 Þ þ 2ðp  zÞð1 þ z02 Þ3=2 r ¼ 0 or more simply:    3=2 rz00 þ z0 1 þ z02 þ 2rp 1 þ z02 ¼0 ð237Þ Here z0 ¼ dz=dr and p- is the pressure in the meniscus measured using a scale of melt column pressure with height for one capillary constant. f0 ¼ p=2  30 . It is important that this equation with the “natural” boundary conditions z0 jz¼0 ¼ tgf0 , rjz¼0 ¼ R0 has a first integral: rz0 ¼ const ð238Þ r2 p þ 1=2 ð1 þ z02 Þ

Fig. 66. Photo of the of the crystallization of a thin rod of LiNb03 (left side) using the EFG method. eqs. (233) and (234) also allows one to calculate of growth angle in the case when crystallization front is not clearly delineated..

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

And accordingly the solution: Z R0 C1  r 2 p pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidr zo ðr0 ; a0 Þ ¼ r 2  ðC1  r 2 pÞ r0

97

ð239Þ

where C1 ¼ Ro p  R0 sinf1 ¼ ro p  r0 sinf0 ¼ const, R0- is the radius of the die(EFG capillary), f1- is a tangent slope angle to the meniscus at its base (die edge). Since it is possible in eq. (233) to use the tangent slope f1 and the radius of any arbitrary point of a meniscus profile curve, so there exists the possibility of determining the growth angle. All linear sizes in (232) are determined using the scale of the capillary constant a and the parameters f1 and ri can be defined by use of a photo of the growth of a thin rod (Fig. 66). The values of f1 and ri can be determined by using some sections of the meniscus profile with its appropriate perpendiculars (the right side of Fig. 66). Since from (231) it follows: p¼

H ri sinfi  rj sinfj ; isj ¼ ri2  rj2 a

ð240Þ

and 30 ¼ p=2  f0 , so: sin f0 ¼

pðr0  ri Þ  rj sin fj r0

ð241Þ

It is possible to use all the linear measurements without taking into account the magnification scale in the photo. In this way the growth angle 30¼0 for Lithium Niobate (LiNbO3). was determined. Also eq. (234) allows one to determine the value of the capillary constant, if one knows the height H of the edge of the die over the melt level in the crucible along with the magnification scale M. The value can easily be determined since the die diameter is well known. Since p ¼ H=a ¼ aMh, where h - is the right-hand side of eq. (234), so after replacement of the dimensionless values ri,j with the corresponding dimensional values it is possible to calculate the capillary constant as a ¼ ðH=MhÞ1=2 . In this way a¼4.35mm for LiNbO3 [60] was determined 8.4. Determination of melt surface tension In practice the determination of the surface tension sm for the melts at high temperatures is especially difficult. The known methods of experimental measurement for sm (or the value of capillary constant): maximum pressure in a bubble, methods with layers or a pendant drop or melt separation from a ring etc., are laborious and require special equipment [61]. However, it is possible to propose a simple and accurate method for the determination of the capillary constant, based on the maximum height achievable before separation of the meniscus [60]. The main idea of this method consists in the experimental determination of the non-dimensional ratio rm =zm , where rm is the radius at the wetting perimeter of any solid probe and zm is the maximum stable height of the meniscus lifted under the free melt surface by this solid probe. The moment of separation occurs at this maximum stable height zm and formally when the second variation of meniscus free energy function d2J changes its sign from positive to negative. The detailed theory of this question is discussed in an article by Zhdanov et al [62]. It is easy to imagine that an ensemble of all the meniscus profile curves can be represented as a field of extremals (stationary points) with an infinitely remote center A(r¼N) in the Cz

98

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 67. Explanation for the determination of the variation task of a stability boundary for the central field of extremals. The height of point A* is the maximal stable height for the extremal.

method (Fig. 67, see also Section 2.1) If an equation z¼z(r,C ) describes the extremals of this field, so it is possible to determine for the parameter C the angular coefficient z0r of the extremal at the central point A. The Jacoby stability requirement for any single extremal z¼z(r,C ) consists in the absence any points on this extremal, which are conjugate to the center A, where it intersects infinitely near to the extremal (Fig. 67). The envelope of extremals is called the C-discrimination line and it determines the stability boundary of extremals (see book p Elsgolz ffiffiffiffiffiffiffiffiffiffiffiffiffiffi[63]). This line is the solution of the Jacoby equation, 0 which looks like h00  rr ½1 þ 3rz 1 þ r 02 h0 þ r12 ð1 þ r02 Þh ¼ 0 for the Cz method. Here the function z¼z(r,C ) is the solution of the Laplace capillarity equation (see above), 2 zðr;CÞ . This equation was very carefully solved numerically and hðrÞ ¼ vzðr; CÞ=vC and h0r ¼ v vCvr one obtained a very detailed values list for the C-discrimination line for Cz menisci in nondimensional coordinates [19]. This unique and single-valued curve zm(rm) (Fig. 68a) is wonderful, since it restricts the maximum stable heights of all kinds of fluids with Cz-like menisci, including high-temperature melts in dimensionless coordinates zm ¼ z=a and rm ¼ r=a. It is possible to represent the tabular data by a single approximation equation [17] Uðr=aÞ ¼ 0:4859ðr=aÞ þ 0:2474 for (r/a>0.2), where Uðr=aÞ ¼ rm =zm and where a ¼ ðsm =rgÞ1=2 is used as the capillarity constant (Fig. 68b).

Fig. 68. (a) e C- line of discrimination in dimensionless coordinates plotted according to data from [19] (b) -graph of the function Uðr=aÞ ¼ rm =zm . The last function is described by equation Uðr=aÞ ¼ 0:4859ðr=aÞ þ 0:2474 for (r/a>0.2).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

99

Fig. 69 shows an experiment to determine the maximum height of the meniscus exactly at the moment of its separation as measured by an electromechanical contact sensor. The relation R=H ¼ U ¼ rm =zm is determined experimentally since the solid probe size R and also the meniscus height H at separation are well known. By using the equation derived above for U(r/a) (or Fig. 68b) and the experimentally determined value of U it is possible to define the value rm (or zm) and then calculate the capillarity constant a ¼ R=rm (or a ¼ H=zm ). The unknown surface tension sm ¼ rl ga2 is calculated for liquids (melts) with known density rl. This parameter rl is determined in the same experiment. It is important to note, that this approach allows one to investigate in the same experiment the dependence of the surface tension sm with temperature and melt composition. This approach also allows one to check the initial melt composition and the equilibrium temperature conditions in a crucible before seeding and the start of crystallization in the Cz and LEC processes. In order to improve the measurement precision it is necessary to take into account the presence of the remaining drop on the end of the probe after separation of the meniscus. This involves the probe edges sharpness and the accuracy of its orientation with respect to the melt surface. The measurement of these heights can be carried out using the drive unit of a standard growth puller to move a solid probe. Of great significance for the modelling and design of feedback control in Cz and LEC processes is the determination of the time constant for a temperature perturbation to relax in the thermal zone. Melt cooling is always a notably slow process following heating after a step wise power change [60].

Fig. 69. In this described scheme there is an additional standard puller arrangement. Experimental scheme set up with an electromechanical contact sensor showing the probe position with respect to the melt level in the crucible and the corresponding signal recording arrangement. 1-control unit of an reversible engine, 2- the engine, 3- the screw, 4- the multiturn potentiometer, 5- the platinum probe, 6-the thermocouple, 7-the compensation voltage source, 8-the amplifier, 9-the recorder, 10- drive unit of a standard growth puller.

100

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

9. Modelling the dynamics of crystallization for a vertical floating zone method (FZ) Here the instability in the dynamics of crystallization look like oscillations of crystal diameter which occur with incremental amplitude (and synchronously with the meniscus height). It is obvious, that there exist conditions of extreme height (and zone length) at which the meniscus will lose mechanical stability and become detached. How can one define this limiting meniscus height for the vertical FZ method? and the next question which is especially important for FZ method - is it possible that loss of mechanical stability is similar to a meniscus satisfying stationary crystallization conditions? The sections below are devoted to the analysis of capillarity related to static and dynamic stability in the FZ method. 9.1. Introduction. Capillarity in the floating zone method The meniscus shape and height in the Floating Zone method is determined as the solution of the Laplace capillary equation which equivalent to the Euler-Lagrange equation for the meniscus freepffiffiffiffiffiffiffiffiffiffiffiffiffi energy function in a gravitational field. This function is ffi Rz F ¼ p 00 ð2sLG r 1 þ r 02 þ gzr 2  pr 2 Þdz for an axially symmetric liquid column with a height z0. The first right-hand member in this integral is the meniscus surface energy, the second item represents the volume part of the general energy and third item is the consequence of the constancy condition for the meniscus volume. The parameter p appears as the indefinite minimum determination for the Lagrange coefficient for the first variation of this function dF¼0 and has a physical sense e it is the pressure in the meniscus. After substitution of the appropriate derivatives in the Euler-Lagrange equation F0r  F00zr0 F00rr r 0  F00r0r r 00 ¼ 0 we receive the Laplace capillary equation r 00 ð1 þ r02 Þ3=2  r 1 ð1 þ r02 Þ1=2 þ sLG ðp  rl gzÞ ¼ 0. This equation in a dimensionless form with the capillary constant a ¼ ð2sLG =rl gÞ1=2 as a unit of linear dimensions is:   3=2  r¼0 ð242Þ r 00 r  1 þ r 02 þ 2ðp  zÞ 1 þ r 02 with - the non-dimensional variables r ¼ r=a, z ¼ z=a and p ¼ p=rl sLG a. If we use the variables z0r ¼ ðrz0 Þ1 and z00 ¼ ðz0 Þ3 rz00 , then the Laplace equation for axially-symmetric menisci in the usual, cylindrical coordinate system Z ¼ ZðrÞ is:    3=2 r¼0 ð243Þ z00 r þ z0 1 þ z02  2ðp  zÞ 1 þ z02 Here we use the simplest notation without the upper line in r; zand p. In eq. (243) the third term is specified with two opposite signs. It is necessary in eq. (243) to correctly assign in the case of convex-concave two-valued menisci with an unambiguous projection of a profile curve on the r axis. The "plus" sign is correct for the meniscus or for part of this meniscus profile curve with a tangent of slope Z0 < 0 and respectively it is necessary to use the "minus" sign for that part of meniscus with Z0 > 0, (Fig. 70). Note, that with coordinates r ¼ r(Z ) such an ambiguity disappears. Different menisci types are considered in detail in reference [64]. For the FZ method pulling in the “upward” or “downward” direction (involving the gravitation vector g) the most interesting case occurs with a cylindrical crystal having a large diameter r0 [a. The Laplace eq. (243) simplifies for axially symmetric menisci at the pulled crystal and having a large radius of radial curvature, in the upward direction:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

101

Fig. 70. Different types of menisci in the FZ method. (a) - Crystallization in an "upward" direction showing the effect of the gravitation vector g. (b)-Crystallization in the "downward" direction. Under stationary conditions r¼r0, the angle to the meniscus and the vertical is equal to the growth angle e0. At the melting front the tangent slope f1 to the meniscus at r¼R0 is variable.

 3=2 z00  2ðp  zÞ 1 þ z02 ¼0 The first integral of this equation is: 1=2  ðp  zÞ2 ¼ C1 1 þ z02

ð244Þ

ð245Þ

Subsequent integration using the boundary condition z(R0)¼0 and z0ðr0 Þ ¼ tg40 gives the elliptical integral which is applicable for both one- and two-valued menisci: Z h  2 ð246Þ C1  ðz  pÞ F1=2 dz rðh; 40 ; pÞ ¼ R0  0

here C1 ¼ cos40 þ ðh  pÞ2 , F ¼ 1  ½C1  ðz  pÞ2  and h-is the limiting height of the menisci at a predetermined value of the growth angle 30. The sub integral function F1=2  0 or 1  cosð40 Þ þ ðh  pÞ2  ðz  pÞ2  1 allows one to determine the permissible and forbidden regions for such menisci. For this case Fig. 71 shows typical calculated plots of the limiting heights and regions for different types of meniscus. From this inequality it follows that the maximal possible meniscus height for pulling upward pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffi is hm ¼ p þ 1  cosð40 Þ with a pressure p ¼ 2 in this meniscus. If this pressure is increased so then the meniscus height is decreased. Similarly for the case of crystallization in the downward direction (to coincide with the gravity vector g) eq. (242) is:  3=2 ¼0 ð247Þ r00 þ 2ðp þ zÞ 1 þ r 02 and for the boundary conditions rðhÞ ¼ r0 , r 0 ðhÞ ¼ ctg40 its solution is: Z h   C2 þ ðz þ pÞ2 F2 dz rðh; 40 ; pÞ ¼ R0  0

ð248Þ

102

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 71. (a) By the use of eq. (246) the menisci heights for the case of pulling upward large-diameter crystals r0 >> a have been calculated. The dotted line corresponds to the forbidden sizes for stationary menisci having a predetermined growth angle of 30¼12 and with a Laplace pressure in the menisci of p¼1.2. (b)- A diagram which shows the permissible and forbidden regions for such menisci. It is possible to determine from the inequality F1=2  0 [64] the boundaries of the permissible and forbidden region sizes for menisci which provide a stationary crystallization process. The capillarity feature in FZ, having an initial state A allows one to retain stability for any perturbation (detail discussed below).

here C2 ¼ cos40  ðh þ pÞ2 , F2 ¼ 1  ½C2 þ ðz þ pÞ2 . The stationary menisci exist in this case if 1=2 2 2 jcosð4 the denominator F2 > 0 in eq. (248). From the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0 Þ  ðh þ pÞ þ ðz þ pÞ j < 1it follows that the p maximum height is hm ¼ jpj þ 1 þ cosð40 Þ for a (positive) pressure in the ffiffiffi meniscus p ¼  2. Since in this case the maximum meniscus height is greater than in the case of pulling upward, so it follows that crystallization in the downward direction is more important in practice for the FZ method. That however is not the only single advantage of crystal pulling in the downward direction. Fig. 72 shows the limiting heights hðDrÞ, where Dr ¼ rðp; 30 Þ  R0 is calculated using eq. (248): Capillarity in the FZ method allows one to pull in the downward direction single crystals of large diameter 2r >> a following their formation by melting poly crystalline rods having a smaller diameter R0 < r. 9.2. Meniscus mechanical stability in the floating zone method This section illustrates the approach to the problem of mechanical stability of static isothermal axially symmetric liquid menisci and is based on work in reference [62]. In general the second variation of the full free energy d2J determines the mechanical stability of the meniscus. The well-known works of Coriell and co-workers [65,66] report a method of p ¼ const finding the conjugated point as zero for the nontrivial solution of the heterogeneous Jacobi equation for a single connected menisci. However, if one uses the general definition of a conjugated point as zero for the first eigenvalue of the second variation of d2J in subspace, selected by the isoperimetric condition (either constancy of the meniscus volume V ¼ const or the pressure in the meniscus) it is possible to find the stability boundary even for multiply connected menisci.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

103

Fig. 72. Limiting menisci heights for the case of pulling large-diameter crystals r0 >> a in the gravity direction (“downward”) calculated using eq. (248) for 30 ¼ 10. Dash-dotted line represents the calculated boundary line for the changeover of the sign of the coefficient Arhc . Stationary crystallization having an initial state A is unstable to any perturbation in the FZ dynamical system (detail in Section 9.3).

Consider, that there exists a liquid zone between two coaxial solid rings having width C1 and C2 in a gravitational field, directed along the z axis in a negative direction. The meniscus length denoted as a and the internal and external profile curves denoted in Fig. 73 are respectively as r1(z) and r2(z).

Fig. 73. (a) Schematic representation the section of a multiply linked liquid zone. (b) Liquid zone (the right-hand of a full liquid zone section) in the form of a tube of large diameter. In both cases the profile curves r1(z) and r2(z)are symmetrical with respect to the axis A. Here as an example is considered the simplest case of multiple linked menisci, typical for the EFG and Stepanov methods.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

104

The function for the free energy J for menisci, represented in Fig. 73a can be written as: Z a h pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi  i J¼p ð249Þ r1 1 þ r1 02 þ r2 1 þ r2 02 þ r22  r12 z dz 0

In order to determine a minimum of the function J it is necessary to take into consideration the condition of the meniscus volume constancy V¼const (a feature of Floating Zone method): Z a  2  V¼p r2  r12 dz ¼ const ð250Þ 0

Note, that in EFG and the Stepanov method menisci have variable volume but with a constant pressure p¼const in meniscus. All linear dimensions in eqs. (249) and (250) use the melt capillary constant as a scale of length. Functions r1(z) and r2(z) minimize the function (249) and take into account the variation restriction condition (9.9). These functions are the solutions of the Laplace equations:  3=2   ri ¼ 0 where i ¼ 1; 2 ð251Þ ri ri00  1  ri 2 þ ð1Þiþ1 2ðz  pÞ 1 þ ri 02 and satisfy the boundary conditions (see Fig. 73): r1 ð0Þ ¼ R1 ; r1 ðaÞ ¼ R2 ; r2 ð0Þ ¼ R1 þ C; r2 ðaÞ ¼ R2 þ C

ð252Þ

Formally the parameter p is represented in (9.10) as the Lagrange factor and has the meaning of pressure in the meniscus. It is possible to define the extremals r1(z) and r2(z) independently. In order for these extremals to yield a minimum in the function J it is a necessary and sufficient that the second variation d2J to be positively definite: Z a d2 J ¼ udz > 0 ð253Þ 0

in space Q2 (will be described bellow). Simultaneously in this task it is necessary to establish the Clebsh condition. In order to further formalize this task let’s introduce L2(0,a) - the space of functions, summed up with the square in an interval [0,a]. Now consider the space, which is a direct product L2(0,a) L2(0,a). Now determine a scalar product hy1,y2i in this space by the formula: Z a ðy1 ; y2 Þdz ð254Þ < y 1 ; y2 i ¼ Here y1 ¼



0

y11 ; y12



  ; y2 ¼ y21 ; y22 ; yKj ˛L2 ð0; aÞ; j; K ¼ 1; 2; ðy1 ; y2 Þ ¼ y11 y21 þ y12 y22

In this space the norm will be noted as kyk ¼ f< y1 ; y2 >g1=2 . Now the space Q2 is defined as a subset in the space L2 L2: 2 Q ¼ fyðzÞ˛C2 ð0; aÞ C2 ð0; aÞ; yð0Þ ¼ yðaÞ ¼ 0; kyk ¼ 1; hy; ri ¼ 0g; where C2(0,a)-is the space of two continuously differentiable functions, r ¼ ðr 1 ; r 2 Þ and r1(z), r2(z)- are the solutions of eq. (251) for the boundary conditions (252). Vector yðzÞ ¼ ðdr1 ; dr2 Þ ¼ ðy1 ; y2 Þ˛Q2 represents the possible variations for the function (9.8) with the isoperimetric condition (9.9). It is possible to write the orthogonality condition as: < y; r >¼ 0

ð255Þ

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

105

This equation means, that the meniscus shape variations do not change the meniscus volume V and is the consequence of constancy condition of the meniscus volume (9.9). If one uses the Schturm-Liouville operator L: d Ly ¼  ðRyÞ þ Py dz

ð256Þ

which acts in the space L2 L2 with the defined domain M ¼ fyðzÞ˛C2 ð0; aÞ C2 ð0; aÞ; yð0Þ ¼ yðaÞ ¼ 0g and by using the matrices R(z) and P(z) in eq. (256) in the form: ! 1 ð1þrr102 Þ3=2 0 1 r2 ð257Þ R¼ 2 0 ð1þr2 02 Þ3=2 ! 00 1 2ðz  pÞ  r102 3=2 0 00 r ð1þr1 Þ ð258Þ P¼ 2 2ðp  zÞ  2 0 02 3=2 ð1þr2 Þ

We obtain the second variation d2J of the menisci free energy functionJ as: d2 J ¼ < Ly; y >

ð259Þ

It follows from eqs. (256)e(258) that the operator L is the direct sum of the SchturmLiouville one-dimensional operators L1 and L2, acting in M according to the formula:   ð260Þ Ly ¼ L1 y1 ; L2 y2 ; y ¼ ðy1 ; y2 Þ˛M Note, that the Clebsch condition, i.e. a positive definiteness of the matrix R(z)>0, is evident and therefore this guarantees the existence of a simple discrete spectrum for the operators L1 and L2. Let lj and 4Kj ðj ¼ 1; .; N; K ¼ 1; 2Þ are the eigenvalues and eigenfunctions of the operators LK respectively. So li(i¼1,.,N) are the set of eigenvalues lKj of the operator L, and for this operator the eigenfunctions 4i - are represented by the vectors ð41j ; 0Þ or ð0; 42j Þ. If mðrÞ is a minimum value of the quadratic functional d2J in the subspace Q2, then in accordance with the minimax theory for eigenvalues the following inequality is correct (see in Lavrentjev’s book [67]): l1  mðrÞ  l2

ð261Þ

If l2<0, so the function d J is certainly not positively definite. Thus the meniscus (mechanical) stability problem is reduced to finding the minimum mðrÞ in Q2. The vectorfunction which derives the minimum function d2 J ¼ < Ly; y > in space Q2, must be an absolute extreme of a function: 2

< Ly; y > mkyk2  2n< y; r >

ð262Þ

For this eq. (262) the Euler-Lagrange equation is: Ly ¼ my þ nr

ð263Þ

The values m- are the eigenvalues of the function d2 J ¼ < Ly; y >in space Q2. It is known from [67,68] that the i-th values of mi(a) reduce with increasing size of a (which is the variable upper limit at integration). The value a*, satisfys one of the equations mi(ai*)¼0 and is the conjugated point to the coordinate origin, i.e. to zero point (Fig. 73). For the non-negativity of

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

106

u (see (253)) in space Q2 it is necessary and sufficient, that there should not exist a conjugated point within the interval [0,a]. Let’s now find the eigenvalues mi(a). The eq. (263) is equivalent to the second kind Fredholm integral equation: Z a yðxÞ ¼ m Gðx; tÞyðtÞdt þ nf ðxÞ ð264Þ 0

with the symmetrical nucleus G(x,t) which represents the matrix Green function of the operator L. N X 4i ðxÞ4j ðxÞ ð265Þ Gðx; tÞ ¼ li i¼1 Ra Here f ðxÞ ¼ 0 Gðx; tÞrdt and for mslithe solution of eq. (264) can be written as: yðxÞ ¼ mn

N X < f ; 4i > i¼1

li  m

4i ðxÞ þ nf ðxÞ

ð266Þ

Since the function f(x) can source wise be represented by the nucleus G(x,t), so it is possible to rewrite eq. (265) as : yðxÞ ¼

N X < r; 4i >

li

i¼1

4i ðxÞ

ð267Þ

The scalar product of the eigenfunctions 4i ðxÞ with eq. (267) gives: < f ; 4i >¼

< r; 4i > li

ð268Þ

After substituting eq. (266) into the orthogonality condition (255) and using eqs. (267) and (268) we derive that m is independent from n and is satisfied in the equation: N X < r; 4i >2 i¼1

li  m

¼0

ð269Þ

This relation (9.27) is the principal relation, sinceR it enables one to calculate the value m for a the known li ; 4i and r. Since d2 J ¼ < Ly; y >  m1 0 ðy; yÞdz and if one knows the eigenvalue spectrum for the operator L and after determination of the first root m1(a) from the eq. (269), it is possible to define whether the meniscus is stable or not stable only by checking the sign of m1(a). The stability is determined by the value m1¼m1(a)0. 9.2.1. Verification of the accuracy of the mathematical approach and a description of experiments with alcohol-water menisci The proposed approach to the stability problem will be illustrated by the determination of the stability boundary hmax for a straight circular liquid cylinder R. It is possible to obtain such menisci only in the state of zero gravity , so the capillarity eq. (251) looks like the following:    3=2 r¼0 ð270Þ rr 00  1 þ r 02 þ 2p 1 þ r 02 Let’s assume in (9.28) r¼r0¼const and we find, that pressure p0and radius of the the meniscus r0 are related as p0¼0.5/r0. For such menisci the classical result hmax/R¼2p is wellknown from literature.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

107

Since such menisci are singly connected, so the relationship (9.27) looks like: N X ðr; 4i Þ2 i¼1

li  m

¼0

ð271Þ

Here ðr; 4i Þ represents an ordinary scalar product in L2(0,a) and r is a profile curve which is the solution of the Euler-Lagrange equation Now let us rewrite eq. (263) as: y00 þ 2

pp0 þ m 2n y¼ pr0 p

ð272Þ

0 ;where The general solution of this equation looks like y ¼ C1 cosðDzÞ þ C2 sinðDzÞ  ppnr0 þm 1=2 D ¼ ½2ðpp0 þ mÞ=pr0  . After the determination the constants C1 and C2 from the boundary conditions y(0)¼y(a)¼0 and from the orthogonality condition, which can be written in this case Ra as 0 y ¼ 0, we obtain the equation:

2  2cosðDaÞ  aDsinðDaÞ ¼ 0

ð273Þ

It is easily seen that m1(a) satisfies the equality Da¼2p or m1 ¼ ð4p3 r0 =a2 Þ  ðp=r0 Þ and from this it follows m1(a)¼0 if a¼2pr0. So we obtain the well known relation a/r0¼2p first determined by Lord Rayleigh [69] for cylindrical menisci in zero gravity. The most prominent and authentic results for the experimental study of the stability problem for liquid columns are reported in work reported in references [70,71]. In these experiments we used immiscible liquids: the oil and alcohol-water solutions of equal density. In [72] the control of the internal pressure in the main meniscus (between pair of dies) was proposed by use of an additional meniscus. However maintenance of the internal pressure in the main meniscus remains an undecided problem as also is the determination of the meniscus surface tension. Fig. 74 shows an advanced scheme which allows one to continuously supervise the Laplace pressure in the complementary (2) and in the basic (1) meniscus. The proposed scheme allows one to vary the curvature 1/R of the lower meniscus (2) and to choose the predetermined pressure in the main meniscus. It allows one to maintain this pressure p¼2sLL/R while increasing the meniscus height by moving the complementary communicating vessel (A). It is important, that in this experiment it is not necessary to determine the capillary constant a¼(2sLL/rlg)1/2 for the alcoholic-water solution in oil since the calculation at extreme heights for all parameters ðr ¼ r=a; z ¼ z=a and p ¼ p=rl gaÞ can be non-dimensional, so that the required parameter rp ¼ r=R is known with great accuracy as the relationship between a pair of linear measurements. For the experimental measurements it is convenient to use the pair pressures which are predetermined by the radius of curvature 1/Ri of the lower meniscus. When Ri¼0, i.e. 1/Ri¼N the lower meniscus has a flat surface, which is parallel to the edges of the lower solid die. The second case e for the greatest possible pressure at which the lower meniscus has the shape of a semi sphere held in place by the edges of the lower solid die. The controlling shape of this subsidiary meniscus is obtained using a longfocus cathetometer. 9.2.2. Comparison of the calculations with experiment It has been established above, that the mechanical stability of the meniscus is determined by the sign of first root of the eq. (269). For the determination of this root it is necessary to solve

108

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

Fig. 74. (a) Scheme for experiments with an alcohol-water meniscus in oil. The liquids have the same density and do not mix. The lower meniscus (2) sets the Laplace pressure in the main meniscus (1). Movement of the complementary communicating vessel (A) by use of drives (C) and (D) allows one to change the height of main meniscus and simultaneously keep the inner Laplace pressure in the main meniscus. (b) Photos of a single meniscus at the moment when it achieves maximum (limit) height and looses mechanical stability (time duration of the breaking-off process w(12)min).

the Schturm-Liouville task for the differential operator L. The coefficients obtained are the functions of the profile curves. These meniscus profile curves (extremals) can be found after integration of eq. (242) using the Runge-Kutt method. In general it is necessary to take into consideration both (meridional and radial) meniscus curvatures and the presence of gravity. In these cases eqs. (242) or (243) do not have simple analytical solutions required for stability analysis or for the numerical calculation of extremals. For such calculations it is convenient to rewrite eq. (243) as the equation set: 8 dz > > < ¼ tg4 dr ð274Þ d4 tg4 2ðp  zÞ > > ¼ þ : dr r cos4 For stationary crystallization the boundary conditions are: the constancy of the growth angle condition 30 ¼ p=2  40 ¼ const, i.e., zjr¼r0 ¼ tg40 or rjz¼z0 ¼ ctgf ; here r0- is the crystal size, and z0 is the meniscus height. This boundary condition must be satisfied with points having an initial r¼R0 and a final r¼r0. The main difficulty in the calculations based on eq. (274) consist in an inadmissible derivative with increasing 40r /Nand correspondingly (intolerable) a decreasing integration step. That difficulty can be overcome with the help of

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

109

changes in the variables by using x ¼ tgð40 ÞlnðrÞ or r ¼ expðx=tgð40 ÞÞ. As a result of this replacement the equation set (274) looks as follows: 8   > tg4 x > > exp < zx ¼  tg40 tg40   tg4 2ðp  zÞ x > >4 ¼  > þ exp : x tg40 cos4tg40 tg40

ð275Þ

The right part of equation for 4x has the order O(1) and at every integration step of eq. (274) it is possible to utilize the simple Runge-Kutt calculation technique. The first step in the calculation consists in the determination of an appropriate profile curve which satisfys the edge condition zr0 ¼ tg40 (for predetermined parameters r0, R,p, 30 and appropriate initial slopes of the profile curves for an initial integration point zR0 ¼ tgf). In the calculations for the twovalued menisci it is more convenient to use the co-ordinates r¼r(z) [64]. In the calculations every meniscus height a is characterized by an appropriate meniscus volume, which can be assumed to be constant in the analysis for meniscus stability at this height. Calculations of the eigenvalues li are based on the Ritz method. This method presumes, that the operator L is positively definite [73]. In the general case, this operator is bottom restricted. This enables one to reduce the problem of finding eigenvalues of the operator L to ~ in the following way. Let searching for the eigenvalues for the positively definite operator L M¼minP(x),(0xa). If we determine the function:  p1 ðxÞ ¼

 PðxÞ þ 1 M0 if PðxÞ  M þ 1 M<0

~ ¼ ðd=dxÞðRy0 Þ þ P1 y is always positive definite. Note, that the eigenthen the operator Ly ~ are the same, and the eigenvalues land ~l are related as: functions for operators L and L (~ l1  M0 l¼ ~ l þ M  1 if M<0

ð276Þ

In the calculation for the system of coordinate functions the following function P was selected pffiffiffiffiffiffiffi ffi ~ operator is 4n ¼ n bðnÞ eK ¼ 2=asinKp=ax. An approximate eigen element 4n of the L K¼1 K eK ðnÞ and an approximate eigenvalue ~lK was determined from the equation ðLn  ~lðnÞ IÞbðnÞ ¼ 0, ðnÞ ðnÞ where bðnÞ ¼ ðb1 .bn Þ and Iis a unit matrix, and the matrix Ln equals to:    ½e1 ; e1 ½e2 ; e1 .½en ; e1      ..........::   Ln ¼    ..........::   ½e1 ; en ½e2 ; en .½en ; en   and   eK ; ej ¼

Z

a 0

~ K ej dx Le

110

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

The approximate eigenvalues of the initial operator L can be found from the relation (9.34). And the approximate values of m(n)are determined from (9.27) for i¼n as roots of the polynomial:          ðnÞ ðnÞ ðnÞ ðnÞ 2 ðnÞ ðnÞ ðnÞ  m  m  m þ C l l . lðnÞ TðmÞ ¼C21 l1  mðnÞ . lðnÞ 1 3 2 n n m     ðnÞ ðnÞ þ . þ C2n l1  mðnÞ . ln1  mðnÞ ¼ 0

Fig. 75. Shows a comparison of calculated data and experimental results (photo on right side) for an alcohol-water solution meniscus within oil (the case of zero gravity). The Laplace pressure in the main meniscus (right) is controlled by the complementary bottom meniscus (scheme on the left side). This pressure remains constant whilst increasing the meniscus height to the point where it has its maximum stable height. The curve of maximum meniscus heights (C- discrimination curve) at constant pressure rop¼const is shown as line P and the thick line V corresponds to the curve of maximum heights at constant meniscus volume. (a)- case of zero pressure in the meniscus rop¼0. Here the dotted line shows for comparison the curve of meniscus heights for a vertical tangent (i.e. equivalent to stationary growth with 3o¼0 ). (b)- the case of maximum Lalace pressure in the meniscus (the shape of bottom meniscus is a hemisphere at rop¼0.58).

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

111

pffiffiffiffiffiffiffiffi P ðnÞ R a where the coefficients Cn ¼ ðr; 4n Þ ¼ 2=a nK¼1 bK 0 rsinðKp=aÞxdx. In the calculations the number of the coordinate functions eK was chosen between 5 and 10. Fig. 75 shows typical calculated profile curves for menisci in a zero gravity state and curves for the highest stable menisci heights. The point 1.0 is the origin of the calculated extremals field. The envelope curve of maximum meniscus heights (that is C- discrimination curve) at constant pressure rop¼const is shown as line P. The thick line V corresponds to the calculated curve of maximum heights at constant meniscus volume. For some the single extremal for the maximum stable height at V¼const is always greater than the maximum stable height at P¼const. On the same graph the menisci heights corresponding to a zero growth angle 3o¼0 are bound by the dashed line (see Fig. 75a). This zero value growth angle corresponds to the maximum menisci heights for stationary EFG crystallization (for the indicated capillarity conditions and 3o¼0 ). In the case of crystal pulling in the upward direction (against gravity Fig. 76a) the comparison of maximum heights for stationary crystallization and for maximum stable heights shows, that the latter are greater. In other words, the dynamic instability in a crystallization process always precedes losses of mechanical menisci stability.

Fig. 76. Calculated menisci “compressed” and “extended” by a gravity field and corresponding to crystal pulling in the “upward”(a) and “downward”(b) situations. In case (a) the dashed line shows the maximum menisci heights for stationary crystallization with a growth angle 3o¼0 and in case (b) the meniscus height for stationary crystallization can attain the stability boundary with 3o>0 . This case is important for the Floating Zone method.

112

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

There is special interest in the comparison heights for stationary crystallization with the calculated maximum stable heights for the case when the vector gravity g coincides with the positive direction of the Z axis, i.e. for crystal pulling in a “downward” direction. Fig. 76b shows the menisci calculated for both cases with the same parameters r0 and p. Note, that in case of “downward” crystallization with any growth angle 3>0 it is possible to reach the boundary of mechanical instability. It is necessary take in to account this feature of capillarity in the Floating Zone crystallization method, since such a “downward” direction is used in practice. The main conclusion - the mechanical instability of a menisci does not prevent the process of stationary crystallization except for the case of “downward” crystallization. In this case it is necessary to calculate the stability boundaries using the meniscus volume constancy V¼const condition. 9.3. Low order model for vertical floating zone method All the features of mathematical modelling and design for feedback control previously considered for the Cz and LEC methods, will be extended to the FZ method. Fig. 77 shows the scheme for FZ. The stationary growth process giving a constant crystal diameter is possible only if the tangent slope to the meniscus at the three-phase line and the vertical lateral crystal surface is equal to the constant growth angle 30(Fig. 77b). Evidently, any variations of this growth regime can change the crystal diameter. In contrast to the Cz and EFG method in the FZ method the Laplace pressure in the meniscus depends on its volume and shape. The meniscus curvature and consequently the Laplace pressure p directly depend on its volume W. This pressure changes by stretching or constricting the meniscus and by adding an additional substance into the meniscus. The Low Order Model for FZ method is based on four laws: i) the mass balance, ii) the growth angle constancy condition, iii) the heat balance at the crystallization front and iv) the heat balance at the melting front.

Fig. 77. Shows a scheme for the Floating Zone method. It is possible to grow the single crystal in the upward (a), or in the downward directions (b) coincident with the line of gravity g. Note, here the origin for the co-ordinates is placed at the maximum temperature T(0) in the melt zone.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

113

In the FZ procedure independently of the method of meniscus heating (inductive, resistive, laser and so on method) there always exists a maximum temperature T(0)>To in the melt zone which is always near to the melting front (Fig. 77) [74,75]. In the FZ method the total meniscus height is h(t)¼hc(t)hm(t) and the crystallization front position hc(t) and position of the melting front hm(t) can change independently. So in general the LOM for the FZ method uses four independent state variables: the crystal radius r(t), the position of the crystallization front height hc(t), the position of the melting front hm(t) and also the meniscus volume W(t). Linearization of conservation laws allows one to write the equation set of first approximation: 8 d_r ¼ Arr dr þ Arhc dhc þ Arhm dhm þ ArW dW > > < _ dhc ¼ Ahc r dr þ Ahc hc dhc þ Ahc hm dhm þ Ahc W dW ð277Þ dh_ m ¼ Ahm r dr þ Ahm hc dhc þ Ahm hm dhm þ Ahm W dW > > : _ dW ¼ AWr dr þ AWhc dhc þ AWhm dhm þ AWW dW We now briefly describe the coefficient Aij in this eq. (277). The thermal conditions for the growing crystal and the part of melt meniscus situated between the crystallization front and the origin of the coordinates (Fig. 77) can be considered similarly to the thermal conditions in the Cz method (see Section 2). The variations of the crystallization front position dhc(t) can be determined from the heating balance at the crystallization front:   ð278Þ ll Gc  ls Gc ¼ rs L V0  h_ c l

s

where V0 is the crystal pulling rate, L is the latent heat of crystallization and Gcl ; Gcs are the temperature gradients in the examined solid and liquid regions. The approach to linearization of this conservation condition and detailed solution of the thermal task were described in Section 2.9.1. The coefficients Ahc r and Ahc hc for the FZ method are equal to the coefficient Ahr and Ahh for the Cz method. It is evident, that the coefficient Ahc hm ¼ 0 since the temperature gradient Gcl does not depend on the position hmand it is possible to take into consideration as a first approximation the independence of dhc from the meniscus volume W and correspondingly Ahc W ¼ 0. The variation of the rate of the crystallization front dh_C(t) is determined from the equation: dh_ C ðtÞ ¼ AhC r dr þ AhC hC dhC þ AhC hm dhm þ AhC T dT þ AhC V dV0

ð279Þ

Note, that the last two items in eq. (279) appear as the result of taking into account of some possible variations of the maximum temperature dT(0)in the meniscus, and of the crystal pulling rate variation dV0 and of the displacement rate of the melting front dVm. It is necessary to consider these variables as external perturbations and also for control in the (T- and V-) channels. The variation of the melting front position is determined from the heating balance condition:   ð280Þ ll Gml  ls Gms ¼ rs L Vm  h_ m here Vm is the rate of displacement of the melting front. Since the diameter of the molten crystal is invariable 2R0¼const, so it is possible to consider that the gradient in the solid phase m Gm s ¼ const. The temperature gradient in the meniscus near the molten crystal Gl is greater, c then Gl in order to support the process of melting of the solid phase and correspondingly

114

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

always hm
ð281Þ

Considering the mass balance condition in the FZ method for case of the crystallization of a cylindrical crystal we obtain the equation for the melt volume change as:      _ ¼ prs R20 Vm  h_ m  r 2 V0  h_ c W rl

ð282Þ

_ ¼ 0 and r 2 V0 ¼ R2 Vm . Since in For stationary crystallization r¼r0¼const, h_ c ¼ h_ m ¼ 0, W 0 0 practice FZ the molten crystal has a cylindrical shape R0¼const, so the perturbations in the system are caused by variations of crystal size r¼r0þdr and by the crystallization h_ c s0 and melting h_ m s0 rates. After linearization of the mass balance condition (9.40) and taking into account the variations dVmand dV0 we obtain: _ ¼ AWr dr þ AW h_ dh_ c þ AW h_ dh_ m þ AWV dV0 þ AWVm dVm dW c m

ð283Þ

Here AWr ¼ 2prs =rl r0 V0 , AWhc ¼ prs =rl r02 , AWhm ¼ prs =rl R20 , AWV ¼ AWhc and AWVm ¼ AWhm . The dependence of the variations of the radius of the growing crystal for stationary crystallization (when b¼0) is determined from the growth angle constancy condition r_ ðtÞ ¼ VC tgbðtÞ (see Section 2). In general, the slope angle of the crystal lateral surface to the vertical is the function b¼b(r,h,W ) (Fig. 77), so the rate of crystal radius variations is :   ð284Þ d_r ¼ Vc b0r dr þ b0hc dhc þ b0hm dhm þ b0W dW Here the angle bðr; h; WÞ ¼ p=2  3  4ðr; h; WÞ near to the crystallization front. Since h¼hchm it is evident, that b0hc ¼ b0hm . It is important to note, that the angle f (Fig. 77) to the meniscus and horizontal at the three-phase line TFL on the melting front is variable and can adjust to various values. From the coefficients in eq. (284) it is possible to determine the solution of the Laplace capillary equation using eqs. (285) and (286). For reasonable deviations of the lateral crystal surface from the vertical one determines these coefficients taking into account, that cos(b)¼1, b0 ¼ 40 and Arr ¼ Vc 40r ; Arhc ¼ Vc 40hc ¼ Vc 40h ; Arhm ¼ Vc 40hm ¼ Vc 40h ; ArW ¼ Vc 40W . The full derivatives dr and dW are:     dr ¼ rh0 40;p dh þ rp0 40;p dp   ð285Þ   dW ¼ Wh0 40;p dr þ Wp0 40;p dp Exclude in eq. (285) the common parameter e the pressure in the meniscus p, and determine the required derivatives:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

115

  Wp0 Wp0 v4 ¼ ¼ 0 0 ¼ D vr h;W r4 Wp  rp0 W40   0 rp v4 40W ¼ ¼ ð286Þ D vW r;h   rp0 Wh0  rh0 Wp0 v4 0 4h ¼ ¼ D vh W;r Since the melt zone volume for both the single- and two-valued menisci in the case of crystal pulling in the upward direction is equal to: Z h 2 ð287Þ Wðh; 40 ; pÞ ¼ pr h þ 2pr zF1 dz 40r

0

Here F1 ¼ ½C1  ðz  pÞ f1  ½C1  ðz  pÞ2 g1=2 , we obtain the private derivatives on the right part of the eq. (286): Z h sin40 F3=2 dz ð288Þ r4 ¼ 2

0

rp0 ¼ 2

Z

h

ðh  zÞF3=2 dz

0

rh0 ¼ ctg40 þ 2

Z

h

ðp  zÞF3=2 dz

0

W40 ¼ 2pr0 hr40 þ 2pr40 Wp0

¼

2pr0 hrp0

ð289Þ

Z

ð290Þ Z

h

h

zF1 dz  2pr0

Z

0

Z

h

þ 2prp

ð291Þ

3=2

ð292Þ

0

zF1 dz  4pr0 0

zsin4F3=2 dz

h

zðh  zÞF2 0

Wp0 ¼ 2pr0 hrh þ pr02 þ 2pr0 hctg40 þ 2prh0

Z

dz Z

h

h

zF1 dz  4pr0 0

zðp  hÞF3=2 dz

ð293Þ

0

It is easy to show, that the denominator D in eq. (286)) is always positive since after substituting in D ¼ r40 Wp0  rp0 W40 of the corresponding derivatives, after several transformations Rh Rh Rh we obtain D ¼ 4psin40 r0 fð 0 z2 F3=2 dzÞð 0 F3=2 dzÞ  ð o zF3=2 dzÞ2 g. The difference in the curly brackets satisfys the Cauchy-Bunyakovskii inequality [41,74], so the inequality D > 0 is always correct. Moreover, since always h>z in this case 40r ¼ r0 W 0 r0 W 0

Wp0 D

r0

> 0 and 40W ¼  Dp < 0. The

determination of the sign of the derivative 40h ¼ p h D h p is necessary for the numerical calculations. In the case of pulling in the upward direction of a crystal with a large diameter r>>a always 4h > 0. Correspondingly, the volume of a single- and two-valued liquid zone in the case of crystallization in the downward direction is equal to: Z h 2 ð294Þ Wðh; 40 ; pÞ ¼ pr h þ 2pr zF3 dz 0

2

1=2 F2 ;

Where F3 ¼ ½C2 þ ðz þ pÞ F2 ¼ f1  ½C2 þ ðz þ pÞ2 g. In this case the derivatives on the right-hand side of eq. (286) are:

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

116

r40 ¼

Z

h

3=2

sin40 F2

ð295Þ

dz

0

rp0

Z ¼2

h

ðh  zÞF3=2 dz 2

0

rh0

Z

ð296Þ

h

3=2

¼ ctg40 þ 2

ðp þ zÞF2 dz Z h Z h 3=2 0 0 0 W4 ¼ 2pr0 hr4 þ 2pr4 zF3 dz  2pr0 zsin40 F2 dz

ð297Þ

0

0

Wp0 ¼ 2pr0 hrp0 þ 2prp0

Z

ð298Þ

0

Z

h

h

3=2

zF3 dz  4pr0

zðh  zÞF2 dz Z h Z h 3=2 0 0 2 0 Wh ¼ 2pr0 hrh þ pr0 þ 2pr0 hctg40 þ 2prh zF3 dz  4pr0 zðp  hÞF2 dz 0

ð299Þ

0

0

ð300Þ

0

As in the previous case one needs to show that the denominator of eq. (286) is always positive and the derivative 40r > 0 and 40W > 0. After the corresponding substitutions the derivative 40h is described as: 8 2 > Z h Z h Z h < 6 2 3=2 4h ¼ D1 pr02 rp0 þ 4pr0 ctg40 ðh  zÞF3=2 dz  8pr ðp þ hÞ z F dz F3=2 dz 4 0 2 2 2 > 0 0 0 : 0 12 39 > Z h = 3=2 A 7 @  zF2 dz 5 > 0 ; The calculation by use of this expression shows that the sign of this derivative 4h can vary in sign in contrast to the previous case. Fig. 78 shows the typical view for the calculated coefficients Aij for the two considered cases. So for the LOM (9.35) for the FZ method it is possible to rewrite it (to a first approximation) as:

Fig. 78. The typical calculated curves of coefficients Aij for the FZ method and for the two considered cases: (1)crystallization in the upward direction and (2)- crystallization in the downward direction.

0

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

117

1 1 0 1 0 10 0 0 0 Arhc Arhm ArW d_r Arr dr B dh_ C B A C B B 0 C Ahc hc 0 0 C hc r B cC B C CB dhc C B Ahc T Ahc V B _ C ¼B C CþB CB @ d hm A @ A @ A @ Ahm T 0 Ahm V A 0 0 Ahm hm 0 dhm _ 0 CAWhC AWVm B þ CAhc r CAhc hc CAhm hm 0 dW dW 0 1 dT B C ð301Þ @ dV0 A dW Here B¼2prs/rlR0V0 and C ¼ prS =rl R20 . The second item of this eq. (301) represents the perturbations and simultaneously the control channel in the FZ method. The dynamic stability of such system is determined by the roots of the characteristic equation detðsI  Ai Þ ¼ 0. According to the Rauth-Hurwitz criterion for this characteristic equation the FZ dynamic system will be stable if it fulfils such an inequality: Ahm hm < 0

ð302Þ

AhC hC þ Arr < 0

ð303Þ

BArW AhC hC > 0

ð304Þ

 BArW Arr  AhC hC Arr ðAhC hC þ Arr Þ þ AhC r ðAhC hC þ Arr ÞðArhC þ CArW Þ 0

ð305Þ

If the liquid zone in FZ is overheated T(0)>T0, so always the coefficients Ahmhm<0, AhChC<0 and the conditions (9.60) and (9.61) are fulfilled and only for an overheated meniscus is the stationary crystallization process possible in the FZ. Since ArW ¼ Vc 40W and always 40W < 0 and B<0, so the inequality (304) is also fulfilled. For the inequality (305) analysis we make use of Fig. 78 where are represented the view of the calculated coefficients for the examined cases (crystallization in the up- and down directions). With certainty the first two items in (305) are positive and if in the third item it fulfils the inequality jArhC j > CArW , so it fulfils the last inequality (305) and as a whole fulfils the RauthHurwitz criterion. Let’s analyse this condition for the case A (Fig. 77a). Using the equations (288)e(293) we determine: CArW ¼ VC D1 pr02 rp0

ð306Þ

 ð307Þ ArhC ¼ VC D pr02 rp þ J þ Y Rh R h 3=2 R h 3=2 3=2 Here J ¼ 4pr0 ctgð4Þ 0 ðh  zÞF1 dz and Y ¼ 8pr0 ðh  pÞ½ 0 z2 F1 dz 0 F1 dz R h 3=2 2 ð 0 zF1 dzÞ . The inequality (305) is fulfilled if the sum is positive J þ Y > 0. For h>p this sum is positive, since the difference in the square brackets satisfies the Cauchy-Bunyakovskii Rh inequality. According to the eq. (290) rh0 ¼ ctg40 þ 2 0 ðp  zÞF3=2 dz < 0. This means, that the menisci of such stable FZ processes belong to the region on the curves of limiting heights, having negative slope rh0 < 0 and the dynamics is unstable if the menisci belong to the region with positive slope rh > 0 (region L-N in Fig. 71a). If the FZ process takes place at h 2ðp  hÞ 0 F3=2 dz, so after replacement ctgð4Þ in Rh Y by the smaller value 2ðp  hÞ 0 F3=2 dz we obtain J þ Y ¼ 2ðp  hÞðhI12  I22 Þ > 0, were  1

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

118

Rh Rh I12 ¼ 0 F3=2 dz and I22 ¼ 0 zF3=2 dz. So again, the stability criterion is fulfilled for menisci belonging to the region only with negative slope. So above it has been shown, that in the FZ process for crystallization in an upward direction a steady growth process is possible for crystals with a smaller diameter then the diameter of a molten crystal r0R0 are difficult since these regions of the limiting heights correspond to the unstable or the forbidden (due to the capillary condition) menisci. But because of the considerable positive difference in the sizes dr¼r0R0>0it is possible obtain crystal growth in a downward direction. In this case the coefficient ArhC has a variable sign (see Fig. 78). If the coefficient ArhC > 0, it brakes the Rauth-Hurwitz rule (9.63) and correspondingly also fails the stability for the FZ crystallization process. In Fig. 72 by the chain line we show the calculated boundary line for the sign change of this coefficient ArhC . Unfortunately, for most menisci which allow one to get large crystal sizes dr¼r0R0>0 fall into the unstable region rh0 > 0 of the limiting meniscus heights. Nevertheless in this case (pulling in the downward direction) the meniscus height is greater than pulling in the upward direction and it is possible to select the sizes of molten crystal for the growth of required crystals with sizes where r0>R0. It is necessary to note, that change of a sign for the capillary coefficient ArhC does not necessarily cause a process of total instability. For the precise determination of this boundary it is necessary to take into account the thermal conditions in (real) FZ processes. 9.3.1. Controllability in the floating zone method Considering the vector U¼(dT,dV0,dVm)T as the control vector, which represents three control channels: (i) when one uses the temperature dT, (ii) the pulling rate dV0 and (iii) the displacement rate of the melting front dVm, it is possible to determine the possibility of controllability in the FZ method having a state vector X¼(dr,dhC,dhm,dW )T (see eq. (301)). After rewriting of LOM in brief matrix notation as X_ ¼ AX þ BU and by use of the classical approach, based on the Kalman controllability criterion, it is possible to show that the controllability matrix Q ¼ ½BjABjA2 BjA3 B does not have the full rank equal to four. This means that the dynamical object is not completely controllable. But in practice for FZ it is possible to achieve the situation where crystallization process takes place for a fixed heating power, i.e. u1¼dT¼0 and for a constant displacement rate of the melting front u3¼dVm¼0. In this remaining case the single control channel u2¼dV0 - the crystal pulling rate and the controlled object can be described by the reduced LOM: 0

1 0 1 10 1 0 d_r Arr dr 0 Arh ArW @ dh_ A ¼ @ Ahr Ahh 0 A@ dh A þ @ 1 AdV0 _ dW AWr þ AWh AWh Ahh 0 dW 1 þ AWh

ð308Þ

In a similar manner we examine the Kalman controllability criterion Q ¼ ½BjABjA2 B, where A and B-are the corresponding matrix of object (9.66) and now determine, that the rank Q¼3 is full. In other words, the FZ process is completely controllable at constant heating power of the liquid zone and for a constant displacement rate of the molten crystal with a constant diameter 2R0¼const. This control is accomplished with a variable crystal pulling rate. The basic features of feedback control design were described above with reference to the Cz and LEC methods and can be directly extended to the FZ method.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

119

The final general observations It seems to the author that non stationarity in Cz and LEC processes is the more significant difficulty compared with nonlinearity. For new crystal diameter control systems it is possible to propose some solutions which should be overcome this problem. But it is the author’s opinion, that the most perspective way will be based on an approach involving modern fuzzy control systems. Acknowledgements Author is grateful to Professor T. Duffar for idea of this review. My especial gratitude to Professor Brian Mullin for the given opportunity to present this material in the journal PCGCM and his steady and friendly attention, also for his huge editorial support. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

V.V. Voronkov, Solid State Physics 5 (2) (1963) 571. V.V. Voronkov, Sov. Phys.-Cristallogr. 19 (1974) 922. V.V. Voronkov, Sov. Phys.-Cristallogr 23 (1978) 249. V.V. Voronkov, Bull. Russian Acad.Sci., Ser. Phys.., 44 (1980) 226 V.V. Voronkov, J. Crystal Growth 52 (1981) 311. V.V. Voronkov, Bull. Russian Acad.Sci., Ser. Phys. 47 (1983) 210. V.V. Voronkov, Bull. Russian Acad.Sci., Ser. Phys. 49 (1985) 2467. A.A. Chernov, Sov. Phys.-Usp. 73 (1963) 277. C. Herring, The phys. of powder metallurgy. McGraw-Hill, NY, 1951, p. 143. C. Herring, Phys. Rev. 82 (1) (1951) 87. D.W. Hoffman, J.W. Chan, Surface Sci. 31 (1972) 368. G.F. Bolling, W.A. Tiller, J Appl. Phys. 31 (1960) 1345. V.V. Voronkov, Sov. Phys.-Cristallogr. 12 (1967) 831. W. Bardsley, F.C. Frank, G.W. Green, D.T.J. Hurle, J. Crystal Growth 23 (1974) 341. E.I. Givargizov, Growth whiskers and lamellar crystals from the vapor. Nauka, Moscow, 1977, pp. 303. G.A. Satunkin, V.A. Tatarchenko, Sov. Phys.-Cristallogr. 30 (1985) 772. G.A. Satunkin, J. Crystal Growth, 255 (2003) 170 S.V. Tsivinskii, Inzh.-Fiz.Zh. Akad. Nauk Belorussk. SSR 5 (1962) 59. K. Mika, W. Uelhoff, J.Crystal Growth 30 (1975) 9e20. D.T.J. Hurle, J.Crystal Growth 63 (1983) 13e17. A. Borshforth, An attempts to test the theory of capillary action, Cambbridge (1883). T.H. Johansen, J. Crystal Growth 84 (1986) 609e620. G.A. Satunkin, J. Crystal Growth 254 (1995) 172e188. V.V. Voronkov, V.M. Pankov, Materials of IX Stepanov Symposium about Shaped Crystal Growth, Leningrad (1986)58e65. I.E. Zino, E.A. Tropp, Assimptotical methods in task thermal conductivity and thermoelasticity. LGU, Leningrad, 1978, pp. 210. E. Kuroda, H. Kozuke, J. Crystal Growth 63 (1983) 278. G.C. Goodwin, S.F. Graebe, M.E. Salgado, Control system design. Prentice Hall, 2001, pp. 900. C.L. Phillips, R.D. Harbor, Feedback control systems, Fourth edition. Printice Hall, Inc, 2000, pp. 609. T. Cariberg, J. Crystal Growth 85 (1985) 32e39. G.H. Hostetter, C.Jn. Savant, R.T. Stefani, Design of feedback control system ISBN 0-03-057593-1, US College edition. SBS College Publishing, NY,, 1982, p. 540. R.G. Lyons, Understanding digital signal processing ISBN 0-13-108989-7. PRENTICE HALL, 2006, pp. 652. W. Bardsley, D.T.J. Hurle, C.C. Joyce, J. Crystal Growth 40 (1977) 13e20.

120 [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75]

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121 W. Bardsley, D.T.J. Hurle, C.C. Joyce, C.C. Wilson, J. Crystal Growth 40 (1977) 21e28. B. Kuo, Digital control system. Holt, ReinhartWinston, Inc, 1980, pp. 446. D.T.J. Hurle, C.C. Joyce, M. Ghassempoory, A.B. Crowley, E.J. Stern, J. Crystal Growth 100 (1990) 11e25. V.V. Denisenko, Modern automation technologies n4. PROSOFT, 2006. http://www.cta.ru pp. 66e74. V.V. Denisenko, Modern automation technologies n1. PROSOFT, 2006. http://www.cta.ru pp. 78e88. V.V. Denisenko, Modern automation technologies n4. PROSOFT, 2006. http://www.cta.ru pp. 86e97. V.V. Denisenko, Modern automation technologies n1. PROSOFT, 2008. http://www.cta.ru pp. 86e99. G.A. Satunkin, S.N. Rossolenko, Crystal Res.Technol. 21 (1986) 1125. I.I. Bronstein, K.A. Semendyaev, Reference book on mathematics. Nauka, Moscow, 1986, pp. 975. R. Izerman, Digital control systems [Russian translation]. Mir, Moscow, 1984, pp. 541. D. Himmelblau, Applied nonlinear programming [Russian translation]. Mir, Moscow, 1975, pp. 536. X. Kvakernaak, R. Sivan, Linear optimal controlling system, [Russian translation]. Mir, Moscow, 1977, pp. 650. R.G. Lyons, Understanding digital signal processing. Pearson Educations, Inc, New Jercey, 2004, p. 652. D. Bendat, A. Pirsol, The applied analysis of casual data, [Russian translation]. Mir, Moscow, 1989, pp. 538. L.R. Rabinder, B. Gould, Theory and application of digital signal processing. Prince-Hall, New Jercey, 1975, pp. 850. V. Cappellini, C.A. Gonstantinides, P. Emiliani, Digital filters and their applications. Academic Press, Inc, 1979, pp. 360. V.S. Gutnikov, Filtering of measured signals, [in Russian]. Energoatomizdat, Leningrad, 1990, pp. 191. K. Brammer, G. Ziffling, Kalman-Bucy Filter, [Russian translation]. Nauka, Moscow, 1982, pp. 198. G.A. Satunkin, A.G. Leonov, J. Crystal Growth 102 (1990) 592e608. G.A. Satunkin, A.G. Leonov, In book process growth of semiconductor crystal and films. Nauka, Novosibirsk, 1988, 122e129. G.A. Satunkin, A.G. Leonov, V.A. Lobatorin, Weight sensor, Bulleten Izobretenii 8 (1988) 127. USSR Patent 1377599 G01 G900. G.A. Satunkin, A.G. Leonov., Weight sensor, USSR Patent 15709609, 22.05, 1989. Ya.V. Vasil’ev, V.N. Mamontov, M.E. Mednik, A.I. Sovertkov, Calculation of weight vertical elastic supprt, Deposit article, CNII TEIPriborostroenija (1980) 22. n1327. V.V. Voronkov, Sov. Phys.-Cristallogr. 19 (1974) 228. A.A. Antonov, V.V. Selin, Electr. Techn. Mat. (in Russian 152 (1981) 20. v3. A. Sanz, J. Crystal Growth 74 (1986) 642. G.A.Satunkin, Proceeding IV International Conference ICSC-2001 (edited by V.P.Ginkin), Obninsk, Russia, Sept.24-28, v2 (2001) 377. G.A. Satunkin, B.S. Red‘kin, V.N. Kurlov, S.N. Rossolenko, V.A. Tatarchenko, Yu.A. Tuflin, Cryst. Res. Technol. 21 (1986) 995. A.W. Adamson, Physical chemistry of surfaces, III edition. John Willey, NY, 1976, pp. 568. A.V. Zhdanov, G.A. Satunkin, R.P. Ponomareva, J. Colloid and Interface Sci. 105 (2) (1984) 334e343. L.E. Elsgolz, Differencial equations and variation calculus. Nauka, Moscow, 1965, pp. 424. G.A. Satunkin, V.A. Tatarchenko, J. Colloid and Interface Sci. 104 (2) (1985) 318e333. S.R. Coriel, M.R. Cordes, J. Crystal Growth 42 (1977) 466. S.R. Coriel, S.C. Hardy, M.R. Cordes, J. Colloid Sci. 60 (1977) 126. M.A. Lavrentjev, L.A. Lusternik, Book in Russian fundumental calculus of variations. vII ONTI, 1935. L. Tclaf, Book in Russian calculus of variations and integral caculus. Fizmatgiz, Moscow, 1966. Lord Rayleigh, Phil. Mag. 34 (1892) 145. J.R. Carruthers, M. Grasso, J. Crystal Growth 13/14 (1972) 611e614. E.A. Boucher, M.J.B. Evans, J. Colloid Sci. 75 (1980) 410e418. M.C. Agafonov, P.A. Gavrilov, L.I. Gubina, V.A. Kislov, V.L. Levtov, Yu.B. Levin, L.V. Leenov, V.V. Savichev, V.A. Tatarchenko, Bull. Russian Acad.Sci., Ser. Phys. 43 (9) (1979) 1935. S.G. Mihlin, Book in Russian “variation methods in mathematical physics”. Nauka, Moscow, 1970. G.A.Satunkin, V.A.Tatarchenko, preprint in russian "Stability Crystallization in Floating Zone Process ", Russian Academy of Science, Solid State Phys. Inst., Chernogolovka 1979, pp. 26. E.A. Brener, G.A. Satunkin, V.A. Tatarchenko, Acta Phys. Acad. Sci. Hungaricae 47 (1e3) (1979) 159e165.

G. Satunkin / Progress in Crystal Growth and Characterization of Materials 56 (2010) 1e121

121

Moscow Institute of Electronic Technique (MIET, Moscow-Zelenograd) 1973. Solid State Physics Russian Academy of Science (Chernogolovka) Candidat phys.-math. Science (Ph.D) 1980. Doctor technical science (Dr.Sci) 1994. 1973e1980 Theory and practice of shaped crystal growth. 1980e1988 Modelling Cz method and design a new sensors and automated pullers. 1989e1990 Invitation work for preperation technology and personnel training of new factory (now “North Crystals”, www.northcrystals.ru). 1991e1998 Modelling and practice of automated LEC crystal growth. 1998e2007 Development the software for automation of plant facilities. 1997e2010 Development of devices and the software for wireless automation, including crystal growing applications.